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ABSTRACT 


This  report  presents  a  combined  analytical  and  experimental  investigation  of 
turbulent  heat  transfer  on  basic  and  composite  configurations  at  hypersonic  speeds. 
The  analytical  results  are  presented  in  Volume  I,  the  experimental  results,  including 
data-theory  comparisons,  are  presented  in  Volume  II,  and  computer  programs  incor¬ 
porating  the  analytical  methods  described  herein  are  presented  in  Volume  III. 

The  two  heat-transfer  prediction  methods  programmed  are  the  laminar-turbulent 
Pr  Mr  momentum  integral  method  and  the  turbulent  nonsimilar  boundary  layer  method. 
This  volume  of  the  report  describes  the  numerical  method  and  presents  flow  charts, 
program  listings,  input  forms,  and  a  description  of  the  output  lor  each  program.  The 
programs  are  written  in  Fortran  IV  language  for  operation  or.  the  IBM  7094. 
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SFCTION  I 


INTRODUCTION 


Two  methods  for  the  calculation  of  turbulent  heating  rates  are  discussed  in 
Volume  I  of  this  report,  the  laminar-turbulent  momentum  integral  method  (PrPr) 
and  the  nonsimilar  method  for  turbulent  boundary  layers  (NSBL). 

This  volume  describes  the  computer  program  for  each  method.  A  description 
of  the  numerical  method  is  given  as  well  as  flow  charts,  program  listings,  input 
preparation,  and  output  description. 

The  first  computer  program  was  written  incorporating  the  Prpr  heat  transfer 
prediction  method  developed  by  R.  A.  Hanks  of  The  Boeing  Company.  The  method 
was  developed  in  the  course  of  the  X-20  program  and  finalized  under  NASA  Research 
Contract  No.  NAS  8-11321  (Reference  1).  In  addition  to  the  PrPr  method  calculations, 
the  program  also  contains  subprograms  for  calculating  pressures,  gas  properties, 
and  streamline  patterns  for  several  special  cases.  The  SRU  1108  version  of  the  pro¬ 
gram  was  written  during  Boeing  Company-funded  research  and  converted  to  the  IBM 
7094  during  this  Air  Force  Flight  Dynamics-funded  study. 

The  program  has  four  major  sections.  Program  A  contains  the  method  per  se. 
Given  the  external  flow  properties  and  wall  condition,  this  section  of  the  program 
computes  laminar  or  turbulent  healing  rates.  All  required  gas  properties  are  com¬ 
puted  by  the  program  from  the  given  pressure  and  edge  velocity.  The  effect  of  tran¬ 
sition  is  estimated  by  matching  laminar  and  turbulent  boundary  layer  momentum 
thicknesses  at  the  transition  point.  The  point  of  transition  is  determined  by  the  pro¬ 
gram  on  the  basis  of  a  transition  Reynolds  number  selected  by  the  user.  Program  A 
can  be  applied  to  any  geometry  for  which  the  external  flow  properties  including 
nonisothermal  wall  effects  can  be  defined  by  the  user. 

The  second  section,  Program  B,  computes  the  local  velocities  and  static  tempera  ¬ 
tures  along  a  streamline  from  a  given  pressure  distribution. 

If  the  streamline  divergence  is  not  known,  Program  C  can  be  used  to  obtain  this 
information.  Program  C  calculates  the  path  of  two  streamlines,  the  streamline  of 
interest  and  one  below,  given  a  two-dimensional  array  of  streamline  angles.  With 
two  streamline  paths  known,  the  divergence  parameter  can  be  computed.  This  calcu¬ 
lation,  however,  assumes  zero  divergence  due  to  lx>dv  geometry. 

The  fourth  section  of  the  computer  program,  called  D,  calculates  or  provides 
information  for  the  three  previous  subprograms  for  four  special  cases;  the  axisym- 
metric  and  two-dimensional  body,  the  hemisphere,  the  swept  cylinder,  and  the  sharp 
delta  wing. 


The  second  computer  program,  the  turbulent  nonsimilar  boundary  layer,  inte¬ 
grates  the  boundary  layer  partial  differential  equations  (described  in  Appendix  C  of 
Volume  I)  using  finite  difference  methods.  The  program  can  calculate  the  turbulent 
boundary  layer  on  unyawed  sharp  or  blunt  axisymmetric  or  two-dimensional  bodies. 
Three-dimensional  flow  effects  are  calculated  using  the  zero  crossflow  pressure 
gradient  approximation  which  implies  no  rotation  of  the  velocity  vectors  within  the 
boundary  layer.  The  program  is  also  limited  to  attached  flow  in  air.  which  is  con¬ 
sidered  as  ideal  gas. 

The  program  is  capable  of  initiating  its  own  boundary  layer  solutions,  given  only- 
external  flow  properties,  for  the  stagnation  point  of  either  sharp  or  blunt  tip  cones 
and  plates. 

Special  inputs  required  are:  pressure,  wall  enthalpy,  a  three-dimensional  flow- 
parameter  and  its  derivative,  the  velocity  gradient  at  the  boundary  layer  edge  and  the 
normal  velocity  at  the  surface,  all  functions  of  the  streamwise  distance  x.  The  user 
must  also  specify  an  initial  and  final  value  of  x,  the  value  of  the  x  increments,  and 
printout  information. 

Output  includes  streamwise  and  normal  velocity,  temperature,  total  and  static 
enthalpy,  shear,  heat  transfer  rate,  a  similarity  parameter,  and  x  derivatives  of 
velocity  and  total  enthalpy.  These  values  are  tabulated  as  functions  of  y  (normal  to 
surface)  at  each  output  position.  Also  printed  in  the  output  are  the  shear  and  heat 
transfer  rate  at  the  wall,  and  the  boundary  layer  displacement  and  momentum 
thickness. 

The  Pr  Pr  progi'am  requires  about  0. 5  to  1  minute  of  computer  time  per  case, 
while  the  nonsimilar  program  requires  2  to  3  minutes  for  a  flat  plate  case.  The  time 
estimate  can  vary  and  depends  largely  on  the  type  of  case. 
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SECTION  II 


THE  LAMINAR-TUIIBULEXT  MOMENTUM 
INTEGRAL  COMPUTER  PROGRAM  (PrHr) 


1.  DISCUSSION 

The  computer  program  desc ribed  in  this  report  is  based  on  the  Hanks  PrPr 
method  of  solution  to  the  boundary  layer  momentum  integral  equation;  and,  as  such, 
contains  all  of  the  assumptions  pertaining  to  the  method  as  discussed  ip  References 
1,  2,  3,  and  Volume  I.  The  program  calculates  aerodynamic  heating  on  arbitrary 
bodies,  with  or  without  streamwise  pressure  gradients  in  laminar  or  turbulent  flow. 
Three-dimensional  flow  effects  are  included  in  the  form  of  streamline  divergence  due 
to  body  geometry  (r)  and  crossflow  or  transverse  pressure  gradients  (f).  The  effect 
of  a  nonisothermal  wall  on  aerodynamic  heating  is  also  calculated  with  modified 
methods  of  Lighthill  and  Seban.  This  is  discussed  in  further  detail  in  Reference  1. 
The  program  described  herein  is  applicable  only  to  air  in  chemical  equilibrium 
(References  4-7). 

a.  Method  of  Solution 


The  PrPr  integral  program  is  actually  four  programs,  any  of  which  can  be  oper¬ 
ated  separately.  The  first  program,  Program  A,  calculates  heating  rates  for  any 
body  shape  at  any  flight  or  ground-facility  test  condition  given  the  external  flow  proper¬ 
ties,  gas  properties  at  the  wall,  gas  properties  at  the  stagnation  condition,  total 
streamline  divergence  and  the  streamline  divergence  due  to  body  geometry  alone. 

Program  A  has  an  additional  capability  of  calculating  transition  effects  based  on 
a  transition  Reynolds  number  selected  by  the  user.  The  program  matches  laminar 
and  turbulent  boundary  layer  momentum  thicknesses  at  the  transition  point. 

In  addition,  the  program  calculates  reference  heat  transfer  coefficients  which  are 
used  to  normalize  the  calculated  local  heat  transfer  coefficients.  The  laminar  heat 
transfer  coefficients  are  all  normalized  In  the  value  of  the  hemisphere  stagnation- 
point  heat  transfer  coefficient  evaluated  at  the  same  free  stream  conditions  as  the 
body  of  interest.  The  heat  transfer  coefficient  on  the  stagnation  line  of  a  f>0°  swept 
infinite  cylinder  is  used  to  normalize  all  turbulent  heat  transfer  coefficients.  The 
hemisphere  and  cylinder  radius  are  items  of  input. 

Program  B  calculates  the  boundary  edge  velocity  and  temperature  given  the 
pressure  distribution  along  the  streamline  or  body. 


Program  C  computes  streamwise  pressure  and  divergence  parameters  given  a 
two-dimensional  field  of  pressures  and  streamline  angles.  Since  pressure  data  and 
streamline  angles  are  often  most  easily  available  as  spanwise  plots  at  various  stations, 
considerable  crossplotting  is  required  to  provide  specific  values  along  the  streamline 
passing  over  any  specified  point  on  the  body.  This  work  is  performed  by  Program  C, 
which  begins  at  the  point  for  which  the  user  desires  to  calculate  the  heating  rate,  and 
traces  out  the  streamline,  proceeding  upstream  to  a  specified  boundary  (e.  g. ,  the 
boundary  layer  origin).  The  streamline  angle  is  determined  by  double  interpolation 
of  a  two-dimensional  array  of  angles  which  is  input.  The  streamline  angle  is  then 
used  to  determine  the  coordinates  of  a  point  on  the  streamline  a  distance  dx  upstream. 
This  procedure  is  repeated  until  the  upstream  boundary  is  reached. 


Tj,  feet 


feet 


It  also  interpolates  as  required  to  obtain  initial  values  at  the  upstream  boundary 
of  the  input  array  (£,  17)  thus  providing  all  information  required  for  Program  B. 

Program  D  provides  the  pressure  and  streamline  arrays  and  initial  values 
required  by  Programs  A,  B,  and  C  for  several  specific  bodies.  These  optional 
configurations  are: 

D-l  Arbitrary  two-dimensional  bodies  at  an  angle  of  attack  and  arbitrary  axi- 
symmetric  bodies  at  zero  angle  of  attack.  The  bodies  may  be  sharp  or 
blunt  and  composed  of  wedges  and  plates  or  cones  and  cylinders. 

D-2  Hemisphere 

D-3  Swept  infinite  cylinder  (0  <  A  <  90°) 

D-4  Sharp  delta  wing  (60°  <  A  <  80°  and  a  >  0) 
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For  options  D-l  and  D-2.  the  user  provides  a  table  of  body  width  or  radius  at 
various  locations  along  the  axis.  For  options  D-2  and  D-3,  the  cylinder  or  hemis¬ 
phere  diameter  must  be  input  and  for  options  D-3  and  D-4,  the  cylinder  or  delta  wing 
sweep  angle  must  be  provided.  For  option  D-4  the  angle  of  attack  must  be  stated. 

In  all  options,  the  user  provides  the  free  stream  conditions. 


Options  D-2  and  D-3  (hemisphere  and  swept  infinite  cylinder)  use  pressure  distri¬ 
butions  stored  within  the  program.  These  distributions  are  shown  in  Figure  1.  Options 
D-l  and  D-4  calculate  pressures  from  the  local  flow  deflection  angle  and  a  modified 
Newtonian  pressure  expression.  For  6  >  0 


2 

sin“  6 


=  1.05  •* 


1.  1025 


(M  sin  6)' 


1/2  2 
1. 27S  sin  6 


M 


0.6 


and  for  6  <  0 


sin"  6 


(1) 


(2) 


This  approximate  analytical  method  was  devised  for  predicting  pressures  from 
the  Newtonian  relationship  in  which  K  is  allowed  to  vary. 


C 


P 


K  sin  6 


(3) 


The  relationship  for  K  was  chosen  to  obtain  the  best  agreement  with  the  exact  solutions 
for  a  wedge,  cone,  and  blunt  body  stagnation  point.  A  plot  of  the  modified  pressure 
coefficients  is  shown  in  Figure  2. 


The  D  options  also  provide  velocity  and  streamline  angle  distributions  along  a 
streamline.  In  option  D-l,  the  velocity  is  computed  from 


and 
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Figure  2:  BOEING  MODIFIED  NEWTONIAN  HYPERSONIC  PRESSURE  COEFFICIENTS 


and 


where  a*  is  the  velocity  at  the  sonic  point 


(6) 


Program  options  D-2  and  D-3  obtain  the  edge  velocity  by  integration  of  the  pres¬ 
sure  gradient  along  the  streamline,  while  D-4  calculates  ue  from 


u 

e 


u 

co 


1  - 


5600 


(7) 


The  development  of  Equation  (7)  is  presented  in  detail  in  Reference  3. 


In  addition,  options  D-3  and  D-4  for  the  swept  cylinder  and  delta  wing  also  provide 
a  two-dimensional  array  of  streamline  angles.  The  streamlines  angles  on  the  swept 
cylinder  are  obtained  from 


tan  0 
e 


u^  sin  A 

OO 


(8) 


where  the  spanwise  velocity  ve  is  obtained  by  integration  of  the  pressure  gradient 
normal  to  the  leading  edge.  Delta  wing  streamline  angles  are  obtained  from  a  corre¬ 
lation  curve  and 


0 

e 


(90  -  A) 


/d 

I0\  1-9 

|9LE-(9°-A»(d 

-)  R 

C  CL 

(9) 


where 


is  the  gradient  of  the  streamline  angle  with 
respect  to  the  ray  angle  e 

The  streamline  angle  distribution  across  the  wing  span  is  a  curve  fit  obtained  by 
matching  the  boundary  conditions  given  by  the  angle  and  the  gradient  at  the  centerline 
and  the  angle  at  the  leading  edge.  The  calculations  of  program  D  are  necessarily  some¬ 
what  approximate  since  exact  flow  field  calculations  are  possible  for  only  a  few  simple 
geometries.  The  approximations  used  in  program  D  are  those  discussed  in  Refex'ence  1 
and  are  simple  yet  reasonably  accurate. 
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Program  D  is  easily  modified  and  extended  as  more  information  becomes  available, 
or  as  certain  specific  shapes  become  important. 

b.  PrPr  Program  Schematic 

The  relation  of  the  four  programs  is  illustrated  on  the  following  page:  however,  it 
should  be  understood  that  all  four  programs  are  coupled  and  function  together  as  a 
single  program  when  desired.  The  symbol  OR  indicates  that  the  required  data  can 
come  from  either  of  the  indicated  sources. 

c.  Nomenclature 

a*  sonic  velocity 

Cp  pressure  coefficient 

f  streamline  divergence  due  to  crossflow  pressure  gradients 

free  stream  Mach  number 

P  local  pressure 

P  stagnation  point  pressure 

o 

P  stagnation  point  or  stagnation  line  pressure 
l2 

r  streamline  divergence  due  to  tody  geometry 

ug  boundary  layer  edge  velocity 

u  free  stream  velocity 

OO  *- 

v  spanwise  velocity  at  the  boundary  layer  edge 

x  distance  along  a  streamline 

6  local  flow  deflection  angle  measured  from  the  free  stream  velocity  vector 

A  total  streamline  divergence,  includes  body  geometry  and  crossflow  pres¬ 

sure  gradients  (rf) 

e  ray  angle  on  a  delta  wing 

6  angular  location 


INPUT 


•  Geometry 

•  Test  condition 

•  Wall  temperature 

•  Transition  Reynolds 

number 

•  Printout  information 


Output  from  Program  B  OR 
input  to  Program  A 

•  Wall  temperature 

•  Pressure 

•  Streamline  angle 

•  Streamline  divergence 

due  to  body  geometry 

•  Streamline  divergence 

due  to  crossflow 
pressure  gradients 

•  Transition  Reynolds 

number 

•  _dge  velocity 

•  Edge  enthalpy 

•  Edge  temperature 

•  Edge  viscosity 


Output  from  Program  A 

•  Heat  transfer  coefficient 

•  Skin-friction  coefficient 

•  Heating  rate 

•  Shear  stress 
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dQ  local  streamline  angle 

A  leading  edge  sweep  angle 

(£,  v)  coordinate  system 

prpr  density-viscosity  product  evaluated  at  a  reference  condition 
2.  Prpr  PROGRAM  DESCRIPTION 

Presented  in  the  following  section  is  a  description  of  the  program  and  subroutines. 
Equations,  nomenclature  input,  and  output  from  each  subroutine  are  presented  in 
alphabetical  order  by  the  subroutine  title. 

Before  proceeding,  a  description  of  the  program  coordinate  system  is  required. 
The  program  operates  in  a  two-dimensional  curvilinear  coordinate  system  where  the 
coordinate  axes  are  designated  as  r?  and  £.  Restrictions  on  this  coordinate 
system  are: 

1.  t)  and  ^  are  orthogonal. 

2.  Both  tj  and  £  must  be  parallel  to  the  surface  over  which  the  boundary  layer 
flow  is  being  considered,  and 

3.  £  is  taken  to  be  in  the  general  direction  of  the  fluid  flow. 

Thus,  no  coordinate  axis  projects  outward  or  normal  from  the  surface  area.  In 
the  case  of  a  flat  surface,  rj  and  £  will  lie  in  the  same  plane. 

CYLINDER 


x,  streamline  distance 


FLAT  PLA'"E 


V 

X 

- 

r 


ii 


a.  Pr^r  Program  Flow  Chart 


TBLP 


TABLE  2 
TABL  IB 
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b.  Subroutines 


The  main  program  of  the  laminar-turbulent  prpr  momentum  integral  computer 
program  —  II  (AS2419)  —  controls  the  input  and  the  logic  flow  to  the  required  mathe¬ 
matical  subroutines.  The  input  logic  is  arranged  to  minimize  the  modification 
required  to  change  the  input  formats  or  the  source  of  inputs  (e.  g. ,  tape  input  or 
direct  coupling  to  another  program). 

1)  ATMOS  (ALT,  SONIC,  PINF,  RHO,  TEMP) 

Given  the  altitude,  ATMOS  calculates  the  free  stream  conditions.  This  subroutine 
is  program  AS1772,  "Standard  Atmosphere  Properties,  1962". 


Input 

ALT 

altitude,  ft 

Output 

SONIC 

a« 

free-stream  velocity  of  sound,  ft/sec 

PINF 

Poo 

free-stream  pressure,  lb/ft2 

RHO 

Poo 

free-stream  density,  slug/ft3 

TEMP 

^oo 

free-stream  temperature,  ®R 

2)  AXI2D 

(Option  D-l) 

AXI2D  calculates  the  flow  conditions  for  arbitrary  axisymmetric  and  two- 
dimensional  geometries  using  the  following  procedure: 

(a)  With  the  geometry  specified  by  the  coordinates  (^j, yj),  the  slopes  and  pres¬ 
sures  are  assumed  to  occur  midway  between  the  input  geometry  points. 


6i+l/2 


tan 
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Where  y  is  the  local  body  coordinate: 


The  end  points  Xj  and  xj  and  the  pressures  at  the  end  points  are  found  by 
linear  extrapolation. 

(b)  If  P  at  either  end  point  is  less  than  P^,  set  the  respective  P  =  P^.  If  P  at 
either  end  point  is  greater  than  PQ,  set  the  respective  P  =  P0. 

(c)  Initialize  x*  =  Xj  and  calculate  the  flow  conditions  at  each  increment  dx  along 
the  streamline 

x*  =  x*  +  dx  i  >  2 
1  i-l 
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Linear  interpolation  is  performed  to  find  P/Pg^fx*)  at  each  x*  using  the 
tables  of  P  and  x  calculated  previously. 

If  P(x,)  <  .  528  P0 


u  (x*)  =  u 
e 


OO 


■7 

I 

) 

t 

i 

I! 


If  P(xj)  >  .528  PQ  find  the  sonic  location  X*  where  P  =  .528  PQ  from 
interpolation  of  the  x  and  P  tables.  Then 


a*x 


(x*)  =  — —  when  x  <  X* 
e  X* 


u  (x*)  =  a*  6  1 

e  I  . 


/  P  \2/? 

IT- ) 

\  o  / 


II 

f 


1/2 


when  x  >  X* 


a*  =  (H/3)1^2  Sonic  velocity 


The  remaining  flow  conditions  are  calculated  from 


A(X*>  =1.0 

i 

—  (x*)  =1.0 
r. 

l 


Two-dimensional  body 


A(X*)  =  y(x*) 
i 

» 

—  (x*)  =  y(x*) 
r. 

l 

i  (x*)  =  H  -  |  u  2(x*) 

e  c  e 

Wg(x*)  =  f^P/PgL(x*),  ig(x*)j 
ZTe(x*)  =  f(p/pgL(x*),  yx*)) 


Axisymmetric  body 


[SOMEGA) 
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(d)  Step  (c)  is  repeated  until  the  flow  properties  have  been  calculated  for  all 
points  on  the  streamline. 

Input 


H 

H 

XSI 

4 

Y 

y(4) 

XI 

xi 

DX 

dx 

POPSL 

VP3L 

ACH 

PINF 

p«, 

total  enthalpy,  ft2/sec2 

coordinate  in  free-stream  flow  directicn,  ft 

geometry  coordinate  orthogonal  to  ft 

streamline  coordinate  at  which  calculations  are 
to  begin,  ft 

streamline  coordinate  increment,  ft 
stagnation  pressure,  atm 
free-stream  Mach  number 
free-stream  pressure,  lb/ft2 


Output 


II 

XSTAR 

POPSLX 

UE 

AIEE 

ZTEX 

OMEGA E 

RORI 

DODI 


x* 

p/pSL(x*> 

ue 

ie(x*> 
(ZT)e(x*) 
we(x*) 
(r/q)  (x*) 

(A/Aj)  (x*) 


number  of  points  calculated  along  streamline 

streamline  coordinate,  ft 

local  pressure  along  streamline,  atm 

edge  velocity,  ft/sec 

edge  enthalpy,  ft2/sec2 

edge  compressibility-temperature  product,  °R 

edge  vircosity  parameter,  slug/ft-sec-®R 

divergence  parameter  due  to  body-shock  layer 
geometry 

distance  between  streamlines,  or  total  streamline 
divergence 
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3)  BEQI 


BEQI  calculates  the  equivalent  distance  parameters  at  the  initial  streamline 
coordinate  Xj.  Following  is  the  sequence  of  operations: 

(a)  Initialize  =  1.0 

x  x 

(b)  Divide  the  first  streamline  increment  dx  into  10  intervals,  dx/10.  Inter¬ 
polate  to  find  the  flow  properties  at  each  interval. 

(c)  Interpolate  for  i  =  f(T  ,  P/P_T)  [IWALL] 

w  w  SL 

ZTS 

=  f  (H,  P/PgL)  [SOMEGA] 

WS 

ZT 

w 

=  f(i,P/P_.)  [SOMEGA] 

cj  w 

w 

(d)  Calculate  the  reference  terms 

w  »  p  u  i  ZT  ,  i  ,  cj  [RORMUR] 
r  r  r  r  r  r 


(e)  Calculate 


d(A/A.) 

dx 


[  DERV] 


(f)  Calculate  ig  and  EXPK  from 


x*  d  (A/A.'; 

_ 

(A/A.)  dx 


If  N  <  .05;  .99  <  N  <  1.01;  EXPK=  0,  ig  SL  =  ie 

-fN(N-l) 

.05  <  N  <  .99;  EXPK  =  -.194e  J  ,  0SL=(N'1)0 

n  /N-l  \ 

N  >  1.01;  EXPK  =  .194  e  >  0SL  "  (n)  6e 
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vP  =  V°s6SL;  *e,  SL  =  H  2" 


(g)  Obtain  E^  and  E^,  (EBAR] 

(h)  Calculate  the  equivalent  distance  parameter 

b  /  b  \  x* 

-f  ,'1FG(F)  xiG(xi)(-fj  *f  G(x)dx 

X*  '  '  X,  X. 


where 


I  I 


gl = 


ana  A/A. 

f/f.  =  — 7-1 
1  r/r. 

1 

The  integration  is  performed  by  applying  Simpson's  rule  to  the  integrands 
Gl  and  G-j-  to  finu  beq}l/x  anc*  ^*eq  t/x>  respectively. 

(i)  bgq/x  is  calculated  for  10  intervals.  The  10th  bgq  /x  is  used  as  the  new 

(b™  /x)„  ,  the  first  interval  is  divided  into  100  intervals,  and  the  above 
eq  xj 

calculations  are  repeated  for  the  first  10  intervals. 

(j )  The  calculations  are  repeated  for  5  iterations  (i)  where 

\m  ]  .  \m  i 

•  10*  l'  x/x,J. 


The  final  (b^  /x)„  is  then  used  as  the  to  /x)v  for  the  general 

l  ew  AII  ii=5  ecl  XI 

calculation  (.QTRANS). 
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Input 

XSTAR 

X* 

streamline  coordinate  at  dx  intervals,  ft 

POPSLX 

p/psL(x*) 

local  pressure  along  streamline,  atm 

TW 

Tw(x*) 

vail  temperature,  #R 

UE 

ue(x*) 

edge  velocity,  ft/sec 

THETAS 

0e(x*) 

streamline  (edge)  angle,  degrees 

OMEGA E 

we(x*> 

edge  viscosity  parameter,  slug/ft-sec-'R 

DODI 

A/Aj(x*) 

dimensionless  distance  between  streamlines 

RORI 

r/rj(x*) 

divergence  parameter  due  to  body-shock 
layer  geometry 

DX 

dx 

streamline  coordinate  increment,  ft 

Output 

BEQXLI 

<beq.L/x)XI 

laminar  equivalent  distance  parameter  at  Xj 

BEQXTI 

<b«..T/x>Xl 

turbulent  equivalent  distance  parameter 
at  xj 

4)  CYLIND  (Option  D-3) 

CYLIND  calculates  the  coordinates  i  and  tj  and  the  two-dimensional  parameters 
0e(|,Tj),  P/Pod,7?)  and  Tw(|,r})  for  an  infinite  swept  cylinder.  The  properties  are 
calculated  as  a  function  of  tj  at  |  -  0  and  are  assumed  to  be  constant  in  the  £  direction. 


Pressure: 


"  ■" 

rpT  I 

P(V  .. 

Em. 

Y 

T2 

P 

0 

Pt 

L  l2  J 

P 

.  o  J 

where 


e  = 


57.3 


V 

rcyl 


P/Prp^(0)  is  found  by  linear  interpolation  on  a  built-in  table  (Figure  1). 
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Stagnation 
Line  Pressure 


M„.  -  M  cos  A 
N  » 


Y'all  temperature  is  found  by  linear  interpolation  on  the  input  array. 
Streamline  angle: 


9  <„  =  tm'1  v(”>  . 
e  u„  sin  A 


where  v(rj)  is  the  tangential  velocity  calculated  in  subroutine  FLOW. 
Then 


P/P^l.7))  =  P/PQ(t?)  for  all  | 

Tw(4.t?)  =  Tw(rj)  for  all  £ 


9  U.V)  =  9  (il)  for  all  £ 

e  e 


(u  )  =  u  sin  A 

e  Xj  00 


Input 

RADIUS 

rcyl 

cylinder  radius,  ft 

SWEEP 

A 

sweep  angle,  degrees 

ACH 

free-stream  Mach  number 

VEL 

free-stream  velocity,  ft/sec 

ETAF 

final  tj  coordinate  for  which  streamline  is  to  be 
calculated,  ft 

XSIF 

final  £  coordinate  for  which  streamline  is  to  be 
calculated,  ft 

TW 

T 

w 

wall  temperature,  °R 
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PINT 

,5co 

<> 

free  stream  pressure,  lb/ft“ 

POPSI. 

*VPSL 

stagnation  pressure,  atm 

Output 

XSI 

K 

coordinate  along  axis  of  cylinder,  ft 

ETA 

coordinate  tangential  to  cylinder,  ft 

M 

number  of  i  values 

N 

number  of  17  values 

ETA  MAX 

^max^) 

maximum  rj  value,  ft 

TII8MX 

^max^  ^ ) 

streamline  angle  at  r)maXt  degrees 

PRESSR 

P/P0«,n) 

pressure  ratio  array 

TWALL 

wall  temperature  array,  *R 

THETAE 

streamline  angle  array,  degrees 

DBTP 

DBTP  is  a 

linear  double 

interpolation  routine.  Search  time  is  reduced 

the  current  search  at  the  search  result  for  the  previous  entry  to  this  routine, 
uses  TBLP  to  perform  the  interpolation. 

Input 

XX  independent  variable  table 

YY  independent  variable  table 

7.7.  dependent  variable  table 

X  independent  variable 

Y  independent  variable 

NX  number  of  values  in  X  table 

NY  number  of  values  in  Y  table 


starting 

DBTP 


NXS  location  in  X  table  at  which  search  is  begun 

NYS  location  in  Y  table  at  which  search  is  begun 

Output 

Z  dependent  variable 

6)  DELTA  (Option  D-4) 

DELTA  calculates  the  body  geometry,  the  streamline  angles,  and  the  two- 
dimensional  pressure  array  P/P0(£,  n)  for  a  sharp  delta  wing.  The  following 
procedure  is  used: 

(a)  If  R  <  0  a  sharp  delta  wing  case  is  indicated  and  6  =  a  (degrees)  for  all  £ . 

(b)  The  pressure  is  calculated  from 

P (i,  *1)  =  |  (-^~)  sin2 6  +  P. 

\sin  6/ 

where 


.  2 

sin  6 


=  1.05  + 


1. 1025  + 


2  2 
sin  6 


1/2 


1.  278  sin2  6 


(c)  The  initial  velocity  is  calculated  from 

<ue>XT  =  uoc|i-  62/5600] 


(d)  The  equations  used  to  calculate  0„(£,  V)  are: 


t  = 


M 


.6 


L  ,  „„  /m 

/tan  <5  V566 

2 

.  -1  1 

2 

l6  + 1-22  (pj) 

(tan  A/ 

+ 

sin  — 
M_ 

oo 

1/2 


Mxt  =  sin  X 
N  * 


1 

90  -  A 


1 

6 


1/2 
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y  -  5($»*)4 

5(***)4+  1 

d*  =  (90  -  A)(l  +  ***)? 

is  found  by  linear  interpolation  on  a  built-in  table  (Figure  3). 


If  ***  <1.2  N  =  {($**) 

v  L 


If  ***  >1.2  N  =  f (A,  ***) 


i  +  ***  _ 
NCL  =  2 -  NCL 

CL  u'f  +  1  CL 

N 


M 


2 

N 


h  -  4  tan  (90  -  A) 
max 


0=  (90-A)NclR+  [0*  -  (90- A)  Ncl|R 
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where 


£  =  tan  (T±/j) 
tan  (90  -  A) 


Then 


8  =8* 
max 


T?  =  TJ.  + 

i  l-l 


(e)  The  delta  wing  option  D-4  assumes  zero  streamline  divergence  due  to  body 
geometry 


Input 

ACH 

M„ 

free  stream  Mach  number 

PINF 

P« 

free  stream  pressure,  lb/ft 

POPSL 

VPSL 

stagnation  pressure,  atm 

VEL 

u„ 

free  stream  velocity,  ft/sec 

ALPHA 

a 

angle  of  attack,  degrees 

SWEEP 

A 

sweep  angle,  degrees 

RADIUS 

R 

leading  edge  radius,  ft 

XSII 

initial  streamline  £  coordinate,  ft 

ETAI 

initial  streamline  rj  coordinate,  ft 

XSIF 

*f 

final  streamline  £  coordinate,  ft 

ETAF 

final  streamline  17  coordinate,  ft 

DELXSI 

A* 

increment  used  to  calculate  £  coordinate,  ft 

DELETA 

At] 

increment  used  to  calculate  rj  coordinate,  ft 

Output 

M 

number  of  £  values 

N 

number  of  rj  values 

XSI 

i 

geometry  coordinate,  ft 

ETA 

■n 

geometry  coordinate,  ft 

ETAMAX 

11  max 

leading  edge  coordinate,  ft 

TH8MX 

e 

max 

streamline  angle  at  leading  edge,  degrees 

THETAE 

*e 

streamline  angle,  degrees 

PRESSR 

P/P0 

pressure  ratio 

URATIO 

u/uw 

velocity  ratio 

7)  DERV  <V,  W,  X,  Y.  Z.  DBYDX,  DX.  ITYPE) 


DERV  calculates  the  derivative  of  a  tabulated  function  f  by  Langrangian  five-point 
formulas.  Given  a  sequence  of  points  x_2»  x_j,  xQ,  Xj,  X2  separated  by  a  common 
increment  h,  the  corresponding  derivative  of  f  at  each  point  is  given  by: 


+  fo> 


y 


y 


Input 

V 
W 
X 

Y 
Z 

ITYPE 
DX  h 

Output 

DBYDX  f’ 
81  DIVERG 


functional  values 


determines  which  derivative  formula  is  to  be  used 
function  increment 


function  derivative 


DIVERG  calculates  P/P0,  Tw,  6 e,  r/r j  and  A/Aj  along  a  given  streamline  by 
the  following  method: 


(a)  Store  the  original  streamline  (1)  data  and  calculate  an  auxiliary  streamline  12) 
such  that 


V 


±  An 


max 
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(b)  Calculate  the  following  at  each  dx  location  on  the  streamline  of  interest 

0  (x*)  linear  interpolation  on  $  (x) 
e  e 

T  (x*>  double  linear  interpolation  on  T^.(  ^ ,  h ) 

P/P  (x*)  double  linear  interpolation  on  P/PJ£,t)) 
o  ° 


-(x*)  =1.0 


/  cos  (**) 

*  n.(x*>  -  n  (x*)  e 

~(x*)  =  — - - - - - 

A.  cos  de  < x *) 

il 

Where  the  subscripts  1,  2,  and  i  indicate  the  original  and  auxiliary  stream¬ 
lines  and  the  point  at  which  calculations  are  begun,  respectively.  If  the 
streamlines  cross  at  any  point,  the  above  calculations  are  begun  a  distance  clx 
downstream  from  the  intersection  point. 


initial  streamline  distance,  ft 


streamline  coordinate  increment,  ft 


FTAF 


final  streamline  coordinate,  ft 


final  streamline  coordinate,  ft 


number  of  £  input  values 


number  of  rj  input  values 


streamline  coordinate  arrav.  ft 


streamline  coordinate  array,  ft 


PRESSR  P/PqU’1!)  pressure  ratio  array 


TWALL 


wall  temperature  array,  'll 


XSIX 


streamline  coordinate,  ft 


ETAX 


streamline  coordinate,  ft 


X 

X 

streamline  coordinate,  ft 

ETA  MAX 

r 1 

max 

n  coordinate  for  geometry  leading  edge,  ft 

THETAX 

streamline  angle,  degrees 

Output 

XSTAR 

X* 

streamline  coordinate,  ft 

PRESRX 

P/PJx*) 

pressure  ratio  along  streamline 

THETAS 

Vx*> 

streamline  angle,  degrees 

DODI 

(A/AjHx*) 

dimensionless  distance  between  streamlines 

RORI 

(r/rjHx*) 

divergence  parameter 

TW 

T 

1  w 

wall  temperature,  *R 

XF 

xf 

final  streamline  x -coordinate,  ft 

II 

number  of  points  on  streamline 

DSIRCH  (Z,  X, 

Y,  ZA,  XA,  YA, 

NZ,  NX,  NY,  IX,  IY) 

DSIRCH  uses  DTAB  to  perform  double  interpolation  on  a  table  of  the  form 

Z  =  f(x,y) 

Input 


X 

independent  variables 

Y 

XA 

table  of  X  values 

YA 

table  of  Y  values 

ZA 

table  of  Z  values 

NX 

number  of  values  in  X  table 

NY 

number  of  values  in  Y  table 

NZ 

number  of  values  in  Z  table 

28 


IX 


order  of  interpolation  in  X  direction 


IV 


order  of  inter|>olation  in  Y  direction 


( iutput 


7. 


dependent  variable 


10)  KBAR  (AYEW,  AYKSL,  H,  SIGR.  EXPK,  POPSLX,  OMEGE,  OMEGS,  GAMC, 
G  A  MO.  EBARL,  EBART) 


EBAR  calculates  the  crossflow  pressure  gradient  parameters  E^  and  E-p  and 
the  pressure  gradient  correlation  parameters  Tc  and  rQ  using  the  following  method 

(a)  Interpolate  to  find 

XT  =  f<P/P__,  1  CT)  [SOM  EGA] 
e,  SL  SL  e,  SI. 

(b)  Calculate  the  local  pressure  in  lb/ft2  and  the  pp  terms 

P 


P  =  P 
r  SL 


(t) 


PTPT  R  Ws 


p  U  ~  |j  U.' 

c  e  R  e 


(c)  Calculate  the  mean  enthalpies 


'avg  2  ^'w  'e.  SL^ 


m,c  avg 
(d)  Interpolate  to  find 


=  i  +  .2  (H  -  i 


e,  SI?  °r 


PTPT 
Pe  ^e 


i  /: 


XT  =  f(P/Pct,  i  )  [SOM EGA] 
m,c  SL  m,c 


ZT  =  f(P/P  i  )  |  SO  MEG  A] 

m.  o  SL  avg 


20 


(e)  Compute: 


(2C  -  .294)  ar 
.4018 


.355 


(I0  -  .294)  <rr 
.4018 


.355 


and 


FKc  ■  («.) 


EXPK 


K 


(2So) 


EXPK 


where 


ZT 


Z  = 


m,  c 


c  ZT 


e,  SL 


ZT 


Z  = 


m,  o 


o  ZT 


e,  SL 


(f )  The  pressure  gradient  correlation  parameters  are  calculated  from: 


E  =  .71764  F 
c 


(g)  Compute  the  crossflow  pressure  gradient  parameters: 


Et  =  1  +  T 
L  c 


/i+  r 

ET  =  (l+.76519  ro)  l—f 
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Input 


A  YEW 

\v 

wall  enthalpy,  ft  /sec 

AY  ESI. 

je,  SL 

t\  n 

edge  enthalpy  at  flow  line  of  symmetry,  ftVsec 

II 

II 

total  enthalpy,  ft^/sec^ 

S1GR 

ar 

Prandtl  number 

POPSLX 

P/PSL 

local  pressure  along  streamline,  atm 

OMEGE 

“e 

edge  viscosity  parameter,  slug/ft-sec-#R 

OMEGS 

“s 

stagnation  viscosity  parameter,  slug/ft-sec- °R 

Output 

GAMC 

rc 

pressure  gradient  correlation  parameter 

GAMO 

rO 

pressure  gradient  correlation  parameter 

EBARL 

PL 

laminar  crossflow  pressure  gradient  parameter 

EBART 

Ej 

turbulent  crossflow  pressure  gradient  parameter 

II)  FIND  (ZA.  Z.  NZ.  XZS) 


FIND  looks  up  Z  in  an  array  ZA.  dimensioned  NZ,  and  stores  the  location  of  Z 
relative  to  ZA(1)  in  NZS.  If  Z  is  not  in  the  range  of  ZA,  the  message  "VARIABLE 
NOT  IN  RANGE  OF  TABLE”  is  printed,  followed  by  the  input  values  of  Z,  Z A ( 1 ) ,  and 
ZA(NZ),  respectively.  NZS  is  then  set  to  1  if  Z  *■  ZA(1)  and  to  NZ  if  Z  ZA(NZ). 
Upon  subsequent  entries  to  FIND  the  search  is  begun  at  ZA(NZS)  where  NZS  is  now 
the  location  of  Z  from  the  previous  entrv  to  FIND.  For  a  normal  exit  NZS  =  i  where 
ZAj  <  Z  <  ZA^j. 

Input 


ZA  array  to  be  searched 


Z 


.‘11 


variable  for  which  location  is  desired 
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Similarly  define 


(7. Tp)  |  as  the  Kth  iteration  for  ZT  at  x*. 

e  x*  I  K  e 


The  iterative  equation  to  be  satisfied  at  each  point  x*  on  the  streamline  is 


|Ux*  Ux*-dx 
K 


<l\R  y  \ 

T  (ZTe)x, 
X 


d_  /_ P 
dx  P 


©I 


K  (P/Po) 

Off 


where 


d(P/P  )  ,  j 

- —  =  -  — (p/pi^  +  -  -  (p/p  » ! 

dx  2  dx(  '  o)  ix*  dx(  W  o' !  x^dx 


(P/Po)  =  2  P~<x*>  +  ^_(x*"dx) 

'  avg  o  o 

The  derivatives  are  calculated  numerically  using  subroutine  DERY. 


For  each  iteration: 


K‘.L 


=  11  -  u  u 


For  the  swept  cylinder  tangential  velocity 


H9  1  ‘ 

‘J  *  H  '  U«T  V  '  2  (u-sin  A)‘ 


Then 


(ZT  c) 
x 

K) 


(x*)  [SOMEGA] 


This  va’ue  of  ZTe  is  then  used  in  the  above  iterative  equation  to  yield  a  new  ux*  . 
These  calculations  are  repeated  at  every  x*  until 


<  € 


where  c  =  10  .  A  maximum  of  ten  iterations  (K  =  10)  is  fixed  by  the  program.  If 

the  relative  error  criteria  has  not  been  satisfied  after  ten  iterations  at  any  point  x*, 
an  appropriate  error  comment  will  be  printed.  If  convergence  at  x*  has  been  estab¬ 
lished  in  the  Ith  step,  then  set  ux*  =  and 


u 


1/2 


I 


I 


A  special  case  results  when  |ue  jx  =  This  corresponds  to  a  stagnation  point 

=  0.  The  iterative  equation  for  ue  becomes 

ue(x*)  =  ue(x*  -  dx)  or  due/dx  =  0.  Because  this  result  is  incorrect,  the  first  step 
is  calculated  using  the  following: 


where  xf  =  0  and  |d/dx(P/P0) 


The  velocity  gradient  at  the  stagnation  point  xj  is  calculated  from: 


where  the  pressure  derivative  is  obtained  from  the  numerical  derivative  equation 


•p_( ~  2  p  (^2,t1^+p 

o  o  '  o 


2 
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The  use  of  this  equation  requires  that  the  first  three  £  input  points  be  equidistant. 


Then 


=  u«  §,  hr)  ix 


(*.)  =H-lK2) 


(ZTe) 


X  ’  \  SL  J 

Ddx  '  X. 


(SOM  EG  A] 


The  calculations  continue,  using  the  previously  described  procedure  where  the 
starting  point  is  now  x}  +  dx. 

Input 

PRESliX  P/PQ(x*)  ratio  of  local  to  stagnation  pressure  along  streamline 

II  number  of  values  in  pressure  array 

U RATIO  (ue/Uoc)  velocity  ratio  at  Xj 

°°  XI 

II  II  total  enthalpy,  ftVsee“ 

VFL  u  free  stream  velocity,  ft/sec 

DX  dx  streamline  coordinate  increment,  ft 

POPSL  P0/P§j  stagnation  pressure,  atm 


Output 


The  following  outputs  are  a  function  of  streamline 
distance: 


edge  velocity,  ft/sec 


POPS  I.X  P/PSL(x*) 


local  pressure  along  streamline,  atm 


AIEE 


ifi  edge  enthalpy,  ft2/sec^ 

ZTEX  ZTe  edge  compressibility-temperature  product,  ®R 

OMEGAE  Wg  edge  viscosity  temperature  ratio,  slug/ft-sec-®R 

XSTAR  x*  streamline  coordinate,  ft 

13)  HEMI  (Option  D-2) 

HEMI  calculates  the  streamline  and  the  pressure  and  divergence  parameters 
along  the  streamline  for  a  hemisphere  geometry. 

(a)  Set  x’J  =  Xj 

(b)  Calculate  x*,  P/PQ(x*),  r/rj(x*),  and  A/A^x*)  for  all  dx  increments  along 
the  streamline  until  0  =  90® 


57.3  x* 

rhemi 


Interpolate  to  find  P/P0(x*)  =  f(0) 


r/r.(x*)  =  sin 


'  x*  \ 

^rhemi  / 


Input 

RHEM  rHEMI  hemisphere  radius,  ft 

XI  Xj  x-value  at  which  calculations  are  to  begin,  ft 

DX  dx  streamline  increment,  ft 

Output 

XSTAR  x*  streamline  coordinate,  ft 

PRESRX  P/PQ(x*)  ratio  of  local  to  stagnation  pressure  along  streamline 
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ROM  r/rj(x*)  divergence  parameter  due  to  body-shock  layer 

geometry 

DOD1  A/Aj(x*)  distance  between  streamlines,  total  streamline 

divergence 

14)  INTIAl.  (LI ,  NTYPE) 

INTIAL  calculates  the  free  stream  conditions  for  wind  tunnel  and  flight  type 
inputs. 

(a)  Wind  Tunnel  Class  i 


=  49 

2 

UO0 

H  =  6000  T  +  — - 

oo  2 

flee  =  0-7  W'J 

=  P./fmeT.) 


(b)  Wind  Tunnel  Class  II 


i_  =  P  - 


when  .*>.90  x  10<,<  i^s  20.3  x  10^  ft“/'sec“  [SONENTJ 


a*  -  .6325  \J  ix  when  i„.<  5.93  x  10^  !l“ ■  sec” 


pw  =  p.AnieT.,) 


P/P  from  fSTACONl 

O  00 


q»  2HT 

oo 

(c)  Flight 

aw  =  f  (altitude) 
=  f  (altitude) 
=  f  (altitude) 


( ATMOS  1 


H  =  6000  T  + 


q#0  =  0.7  P„mJ 


For  wind  tunnel  I  and  flight 

VP«  =  f<M*.  P„.  UJ  [STACON] 


VPSL  -  1 

/P  \  P„ 

I  o  \  * 

\  P_  /  PC.T 
'  */  SL 

Input 

Wind  Tunnel  Class  I,  NTYPE  =  WT1 

ACH 

M 

OO 

free  stream  Mach  number 

TEMP 

T 

x00 

free  stream  temperature,  ”R 

PINF 

p_ 

oo 

free  stream  pressure,  lb/ft^ 

Wind  Tunnel  Class  II.  N’T Y PC  -  WT2 


VEL  u  free  stream  velocity,  ft/sec 

2  o 

II  II  total  enthalpy,  ft  /sec 

POPSL  PQ/Pgj  stagnation  pressure,  atm 

PINT  Pw  free  stream  pressure,  lb/ft" 

Flight,  NTYPE  =  FLIGHT 

ALT  flight  altitude,  ff 

VEL  uw  free  stream  velocity,  ft/sec 

Output 

ACM  M*  free  stream  Mach  number 

H  II  total  enthalpy,  ft2 /sec2 

PINT  P^  free  stream  pressure,  lb/ft 

POPIXF  P0/P*j  stagnation  to  free  stream  pressure  ratio 

POPSL  stagnation  pressure,  atm 

QIXF  qK  free  stream  dynamic  pressure,  lb/ft2 

SOX1C  speed  of  sound,  ft/sec 

VEL  u^  free  stream  velocity,  ft/sec 

15)  IWALL  (TW,  POPSLX,  AYEW) 

1WALL  calculates  the  enthalpy  at  the  wall  as  a  function  of  the  wall  temperature 
and  the  local  pressure. 


(a)  Set  =  Tw/l,8  (conversion  to  6K) 

A 


(b)  Calculate  wall  enthalpy 


i  -  ,  432  T  when  T  <  ::003K 
w  "K  \ 


i  =■  .432  T,  -•  3. 32  x  10  ;’/T  .  -  300\2  when  300  <  T  <  l’iOO’K 

w  "iv  1  "K  '  WK 

iw  =  f (log.  0  P/Pgl,  Tw  );  [TABLE  g|  when  T^  >  18003K 


(c)  Convert  i  v  to  ft2/sec2 


i  I-"-]  -  23031  i  (8jS) 
w  1  2  I  w  \  lb  / 

\  sec  / 


Input 

TW 

Tw 

wall  temperature,  °R 

POPSLX 

P/pSL(x*) 

local  pressure  along  streamline,  atm 

Output 

A  YEW 

*w 

wall  enthalpy,  ft2/sec2 

16)  JAYELL  (A YEW,  H,  AYEE,  ZTE,  ZTS,  S1GR,  BETAS.  POPSLX,  XJL) 

JAYELL  calculates  the  laminar  pressure  gradient  effect  correlation  parameter  JL 
from  the  following: 


(a)  Calculate  the  mean  enthalpy 


i  =  -  (i  +  H) 
m,  s  2  w 


(b)  Interpolate  to  find: 


ZT  =f(P/P_.,i  ) 

m,  s  SL  m,  s 


[SOM  EG  A] 


(c)  Calculate 


(I  -.294) 
s 

.4018 


(¥) 


.355 
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where 


Input 

A  YEW 

\v 

.  2 

wall  enthalpy,  ft'Vsec 

H 

H 

2  2 

total  enthalpy,  ft  /sec 

AYEE 

*e 

edge  enthalpy,  ft  /sec"1 

ZTE 

ZTe 

edge  compressibility-temperature  product, 

ZTS 

ZT 

s 

stagnation  compressibility-temperature 
product,  °R 
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SIGR 


Prandtl  number 


a 

r 

BETAS  streamwise  pressure  gradient  parameter 

POPSLX  P./Pg[  (x*)  local  pressure  along  streamline,  atm 
Output 

XJL  laminar  pressure  gradient  effect  correlation 

parameter 

17)  MUZERO  (OMEGR,  ZTR,  AYER,  AYEE,  H,  POPSLX.  XMUO,  XL) 

MUZERO  computes  the  reference  stagnation  viscosity  p0  and  the  diffusion  effect 
parameter  Z . 

(a)  Interpolate  to  find: 

-p=  f(log10  P/PSL,  ie)  [TABLE  7] 
e 


If 


<  °  ’  “•  V1,  "  # 


(b)  Calculate  and  Z 


H  =  w  (ZT  )  — 
o  r  r  1  ' 


(if2 


ZT 


200 


L(i)zv 


200 


Z  = 


i 

e 


H 


Input 

OMEGR 

reference  viscosity-temperature  r'atio,  slug/ 
ft-sec-  ®R 

ZTR 

ZTr 

reference  compressibility-temperature 
product,  8R 

AYER 

*r 

reference  enthalpy,  ft“/sec^ 
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AYEE 

i 

e 

edge  enthalpy,  ft2/sec2 

H 

H 

total  enthalpy,  ft2/sec2 

POPS  LX 

P/PSL(x*) 

local  pressure  along  streamline,  atm 

Output 

XMVO 

Mo 

reference  stagnation  viscosity,  slug/ft-sec 

XL 

£ 

parameter  indicating  diffusion  effect  on  heat 

transfer 


18)  QTRANS  (Program  A) 

QTRANS  calculates  heat  transfer  data  for  any  body  shape  given  the  external  flow 
properties,  gas  properties  at  the  wall,  gas  properties  at  the  stagnation  conditions, 
and  the  streamline  divergence  parameters.  The  sequence  of  calculations  is  as  follows 
(names  in  brackets  are  subroutines): 

(a)  Set  up  x-printout  array 

(b)  Initialize  /3g  =  0.5 

(c)  Interpolate  for 


V  =  f(Tw'  P/PSL>  llWALLl 


(d)  Interpolate  to  find 


ZT 

s 


Cl) 


S 


f(H,  P/PgL)  [SOM EGA] 


ZT 


w 


W 


W 


■  (<V  P/PSL) 


[SOM  EGA] 


(e)  Calculate 


u; 


P  P  ’ 
rr  Kr 


ZT 


[RORMUR] 
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(f)  Calculate 


du  d(A/A.) 
e  l 


dx  ’  dx 


(DERV] 


(g)  Calculate  i  _T  and  EXPK  from 
e,  SL 


N  = 


*  d(A/A.) 
x*  _ _ r_ 

A  dx 


If  N  5  .  05  ;  .  99  <  N  <  1.  01  ;  EXPK  =  0 ,  i  OT  =  i 

e,  SL  e 

4n(N-1)  _ 

.  05  <  N  <  . 99 ;  EXPK  =  194e  ,  fl  =  (N  -  1)  9 

SL  e 

~3  -  /X  -  1  \ 

N  >  1.01;  EXPK  =  .  194  e  J  .  t>ar  =  ---  g 

SL  \  N  /  e 


(h)  Obtain  and  E-j.  [EBAR1 

(i )  Calculate  the  equivalent  distance  parameters  defined  by: 

b 


eq  _  i 
x  x*  G(x*) 


/b  \  x* 

xiG(xi)  Hr/  +  /  G(x)dx 

'  'X.  *x. 


1 

x*G(x*) 


(b  \ 

c.  G(x.)  +  f 

I  X  \  X  /  J 

V  •'v 


x*-2dx 


X.  X. 


G(x)dx  + 


f 

*v*- 


x*-2dx 


where  the  integrand  G  is  calculated  from 


/  \2  / ,  \2Et 

GL  ■  (r)  (t) 


/  r  \5/4  /  f  \<5/4>ET 
Gt  =  (PrPr)ue /— -j 


and 


f.  r/r. 

l  l 


(x)dx 


44 


■45 


V 


and 


(k)  Calculate  the  reference  Reynolds  number  R  •  calculate  J  [JAVELL] 
and  u  £  (MUZERO]  r’^ 


then 


(VX)L 


•  F 
2  ’  x,Q 


(beq/x). 


;  R 


x  /  .9/16  ’  r,Q 

1. 


/b  \  11/3 
17/20^begj 


L  J 


p  u  u  x  /x  )  x* 
r  *r  e  v  eq  / 


2  2 

M  F 

o  x,  Q 


/*  \  (beq/X). 

l^L  ■  ' 


4/10 


if  Rr,q  a  ^transition  <input)  a  new  value  of  (beq/x  )  is  calcuiated  using  the 
following  equations: 


R 


L,  S 


R, 


p  u  u  S 
rr  rr  e  eq,  L 


p  m  u  S 
_  Mr  r  e  eq,T 


T,  S 


4C> 


R - 

e,L  (x 

e 


,664 


(\s) 


1/2 


pu  9, 
e  e  L 


R 


V2135 


P«9t 

eel 


9,  T  /,  _  >2. 64  p 

"eKoVs*-407) 


'transition 


./Vk.r'4  1/ 

VW .  Jo 


transition 


L  transition 


GTdx 


/b  \  /b  \  C_D 

eq  J  eq  TR 


l- 


x  /T  \  x  /T  GTX* 


If  Rr  Q  <  ^transition  (inPut)»  a  check  is  made  at  this  point  for  x-printout 
values.  If  a  printout  point  does  not  lie  in  the  open  interval  (x*  -  2dx,  x*  +  2dx), 
where  x*  is  the  current  value  of  x,  the  sequence  (c)  through  (k)  is  repeated 
until  one  is  encountered. 


(1)  When  a  printout  point  lies  in  the  above  interval  and  when  the  transition  point 
is  found  the  desired  heat  transfer  data  is  calculated  using  the  following 
equations: 

Calculate  JL  and  nQ  X  [JAYELL]  and  [MUZERO] 


Then 
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Reference  Reynolds  numbers: 

PH  u  (x  /x\  x* 
r  r  el  eq  /. 


R 


r,Q 


2  2 

P  F  ^ 
o  x,Q 


p  |i  u  (S  /x )  x 
r  r  e  \  eq  L 


R 


r.S 


2  2 

u  F 
o  x,  S 


TT  _  MoFx,  SJL 


3/2 


/S  /x  \  > 

(  ^  jL 

3/2 


185  R 


,S 


loS10(Rr,S+3000) 


=  .332 


_ %  Fx,  s(Rr,s) 


1/2 


Skin  friction  coefficients; 


CF  T 
t,c,L 


u  p  u 
e  e  e 


C 


f.e.T 


1 


p  u 
e  e 


where 


P 

e 


P/P 


1.232 


SL 


ZT 


.534 


48 


Heat  transfer  coefficients  based  on  enthalpy; 


-  3/10  „  /  \l/2 

.  332  g^J  p  F  (R 

H  _  _ L  o  x,Q<  r,  Q/ 

L  .  645  /x  /x  \  x* 


„  i>  T  3/10 
gjCjL  MoFx,Q 


.  185  R 


T  .645 


°r  (xeq/x)_x*  |1Ogl0(Rr,Q+  3000>| 


2.584 


hl=  h- 


H„,  =  — 
T  H 


ref,  T 


Nonisothermal  wall  parameters; 


*local  ( *w. 


V  f  A  V- 

ml  l  w  V  /  K 


,  .  W  /  ~  \  w.  w.  ,  /  S. 
local  o/  i=l  \  i  1-1/  i 


where 


_  f,  /§*',3/4 

a  -  1 1  (SL/. 


laminar 


v  -rcv 


turbulent 


\  ~  —  / »T  -  ST  v 
L.  3  L.  L.  -x. 

i  i  i-l/  „  _  f  i 

CCT-’  Li  ’  /° 


G  dx 

ij 


PC- 'Vt  ' 
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Adiabatic  wall  enthalpy: 


i  =  i  a  1/2  (II  -  i  ) 
aw,  I.  e  r  o  e 


i  -  i  +  a  * 1  ^  (II  -  i  ) 
aw.  Ter  o  e 


Nonisothermal  heat  transfer  rates: 


H  /i  -  i 
L  I  aw,  L 


+  ♦ 

w  L 


local 


j'aw,  T  ~  \v  *T 


local 


) 

) 


Isothermal  heat  transfer  rates: 

q  .  -  II,  /i  -  i  \ 

L,  iso  1,  \  aw,  L  w / 


^T,  iso  ('aw, T  *  'w) 


(m)  Print  out  iieat  transfer  data. 


In  order  for  the  sequence  of  calculations  to  reach  this  point,  at  least  one 
printout  point  Xp()  must  lie  in  the  interval  (x*  -  2dx),  (x*  +  2cL\)  where 
x*  -  Xj  ■+  (2k)dx  (INTKGFK  k  ?  0)  is  the  current  value  of  x.  Answers  of 
interest  are  now  printed  using  linear  interpolation  for  those  printout  points 
satisfying  x*  -  2  dr  <  Xp(y  x*.  If  a  printout  point  satisfies  the  criterion 
x*  <  Apo  <  fx*  +  2.ixl  then  all  quanti  ties  are  saved  for  the  current  value 
of  x  to  use  for  the  interpolation  scheme  at  x  +  dx.  Steps  (a)  through  (ml 
are  repeated  until  all  printout  points  have  been  exhausted. 


Input 


Streamline  Functions 


XSTAR  x*  streamline  coordinate  at  dx  intervals,  ft 

POPS  LX  P/Pgf(x*)  local  pressure  along  streamline,  atm 

TW  Tw(x*)  wall  temperature^ 0  K 

A  IFF  ie(x*)  edge  enthalpy,  ft^/sec^ 


GO 


THETAS 

*e<X*> 

streamline  edge  angle,  radians 

DODI 

A/Aj(x*) 

divergence  parameters 

RORI 

r/r^x*) 

UE 

ue<x*> 

edge  velocity,  ft/sec 

ZTEX 

ZTe(x*) 

edge  compressibility-temperature 
product,  “R 

OMEGAE 

we<x*> 

edge  viscosity  temperature  ratio,  slug/ft-sec- 

Initial  and 

Reference  Conditions 

H 

H 

total  enthalpy,  ft2 /sec 2 

XI 

XI 

initial  streamline  value,  ft 

RTRANS 

^transition 

transition  Reynolds  number 

SIGR 

or 

partial  Prandtl  evaluated  at  ZTr 

DX 

dx 

streamline  coordinate  increment,  ft 

XPO 

xPO 

printout  array,  ft 

KF 

KPO 

number  of  values  in  printout  array 

AYEvVO 

i\vn 

reference  wall  enthalpy,  ft2/sec^  (needed 

only  if  Tw  ^  constant) 

BEQXLI 

/*>  a\ 

initial  value  of  bgq/x  (laminar  and  turbulent) 

BEQXTI 

J 

Output 

\  X/x, 

XSTAR 

X* 

streamline  coordinate,  ft 

POPSLX 

P/PSL(X*) 

local  pressure  along  streamline,  atm 

UE 

ue(x*) 

edge  velocity,  ft/sec 

A/Aj 

A/Aj(x*) 

distance  between  streamlines 

ROM 

r/rj(x*> 

divergence  parameter  due  to  body-shock  layer 
geometry 

TW 

Vx*> 

wall  temperature,  °R 

ZTKX 

ZTe(x*) 

edge  compressibility-temperature  product,  °R 

OMEGAE 

we(x*) 

edge  viscosity  temperature  ratio,  slug/ft-sec-0 

AIEE 

edge  enthalpy,  ft2 /sec 2 

A  YEW 

iw(x*) 

wall  enthalpy,  ft2/sec2 

BEQXL 

beq,  L/x<x*> 

ratio  of  laminar  equivalent  distance  parameter 
to  distance  x  along  streamline 

BEQXT 

beq,  T/x(x^ 

ratio  of  turbulent  equivalent  distance  parameter 
to  distance  x  along  streamline 

AYEAWL 

‘aw,  l/x*) 

laminar  adiabatic  wall  enthalpy,  ft2/sec2 

AYEAWT 

‘aw,T<x*> 

2  2 

turbulent  adiabatic  wall  enthalpy,  ft  /sec 

HL 

Hl(x*)/Ho 

laminar  heat  transfer  coefficient  ratio, 
lbm/ft2/sec 

HT 

HT(x*)/Href  T 

turbulent  heat  transfer  coefficient  ratio, 
lbm/ft2/sec 

QLISO 

^iso,  L<x*> 

laminar  isothermal  heat  transfer  rate, 
Btu/ft2-sec 

QTISO 

<liso,T<x*) 

turbulent  isothermal  heat  transfer  rate, 
Btu/ft2-sec 

QDOTL 

qL(x*) 

laminar  heat  transfer  rate,  Btu/ft2-sec 

QDOTT 

qT(x*) 

turbulent  heat  transfer  rate,  Btu/ft2-sec 

CFEL 

Cf,e.L<x*> 

laminar  skin  friction  coefficient 

CFET 

Cf,c,T(X+) 

turbulent  skin  friction  coefficient 

RRQ 

Rr,Q 

heat  transfer  reference  Reynolds  number 

RRS 

Rr,  S 

skin  friction  reference  Reynolds  number 

19)  REF 


REF  i.  i.culates  the  stagnation  heat  transfer  coefficient  based  on  a  reference 
hemisphere  of  radius  RQ;  and  the  reference  heat  transfer  coefficient  based  on  a  60° 
infinite  swept  cylinder  given  the  reference  cylinder  radius.  The  following  procedure 
is  used  to  calculate  the  stagnation  conditions: 

(a)  Initialize  ig  =  H 

(b)  Obtain  the  following  from  the  indicated  subroutines: 


H> 

[IWALL) 

ZT  ,  ZT  ,  ZT  ,  v  ,  u  ,  u> 
w  e  s  w  e  s 

[SOMEGA] 

w  .  .  ZT  ,  i  ,  a 

r  r  r  r  r  r 

(RORMUR] 

JL 

[JAYELL] 

tV  1 

[MUZERO] 

(c)  Calculate  the  stagnation  heat  transfer  coefficient  and  heat  transfer  rate: 


H 

o 


.  664  g 

.645 

a 

r 


R 


HE  MI,  ref 


The  reference  conditions  are  calculated  from  the  following: 

(d)  Reinitialize  M^t  =  1/2  Mw 


(e)  Obtain  the  following  from  the  indicated  subroutine: 

PST/P»'  W  "mN  ISTAC0Nl 

Hv,  ref  |IWALL' 

ZV  ZTe’  ZTs'  V  “V  “s 


(SOMKGA  ] 


(RORMUR] 


ll?  ,  p  p  ,  ZT  .  i  ,  ct 
r  r  r  r  r  r 


r  ,  r 

c  o 


p  .  £ 
o 


(KBAR] 
(MUZ FRO] 


(f)  Calculate  the  following: 


.866  R^...  , 

_ C\  L,  ref 

«1-L  (,+  rc)uMN 


/  1  +  Fo 

Fx,Q  "  1,6  \1  +  .7625  T 


i+  r 


i  +  r 


.  866  p  p  u  x 

,  r  r  *  eq,L 

r.Q.ref  2  2 

“o  FX,  Q 


,332gj£  Mo  x,Q 


ref,  L  .  645  x 

a  eq,  L 


/R  V 

x  ,  (r.Q.ref/ 


n 

.  185  g  £  ^ o  x,  Q 


r,Q,  ref 


ref‘T  ar  645  ^  |IOg10(rtr,Q,ref+30P0) 


2,584 


„  =  H  -  /l 

-  a  U2  \ 

i  (H  -  i  ) 

aw,  ret,  L  \ 

r  f 

,  ^  =  H  -  /l 

1/3  v 

-  cr  ! 

(H  -  i  ) 

aw,  ref,  T  \ 

r  ) 

^ref,  L  * * ref,  L 

(i  T  - 
aw,  L 

\v,  ref^ 

n  —  h 

qref,T  ref,  T 

(law,  T  " 

i  ,) 

w,  ref 

1 


fi45 


0„„  *645 

•  866  u  (T 

oo  j- 

H 

ref,  L 

g£ 

ref, 

*645 

.  8b6  u  a 

OO  j» 

Tref,  T 

13  _  _ 

ref,  T 

Tref 

f,  c,  ref 

Input 

PINF 

Poc 

free  stream  pressure,  lb/ ft^ 

ACH 

free  stream  Mach  number 

VEL 

UO0 

free  stream  velocity,  ft/sec 

VELP 

» 

UMN 

modified  Newtonian  velocity  gradient,  l/sec 

H 

H 

total  enthalpy,  ft2/sec2 

POPSL 

Po/PSL 

stagnation  pressure,  atm 

TWRFF 

T\v,  ref 

reference  wall  temperature,  *R 

RHMREF 

RHEMI.  ref 

reference  hemisphere  radius,  ft 

RCLREF 

R 

CYL,  ref 

reference  cylinder  radius,  ft 

Output 

AYEWO 

»w 

o 

stagnation  wall  enthalpy,  ft2/sec2 

HO 

Ho 

stagnation  heat  transfer  coefficient, 

lbm/ft2-sec 

QDOTO 

% 

stagnation  heat  transfer  rate,  Btu/ft2-sec 

HREFL 

11  ref,  L 

laminar  reference  heat  transfer  coefficient, 
lbm/ft2-sec 

HREFT 

Kref,  T 

turbulent  reference  heat  transfer 

....  ..  •> 

coefficient,  sec 


A  YAW  It  I. 

'aw,  ref,  1. 

laminar  reference  adiaoatic  wall  enthalpy, 
It“/sec“ 

AYAWRT 

'aw,  ref,  T 

turbulent  reference  adiabatic  wall  enthalpy, 
ft2/sec2 

TAl'RI. 

Tref,  L 

laminar  reference  shear  stress,  lb/fi 

TAIJRT 

Tref,  T 

Q 

turl)ulent  reference  shear  stress,  lb/ftr 

QRKFL 

f*ref,  L 

laminar  reference  heat  transfer  rate, 
lbm/ft2-sec 

QREFT 

^rel,  T 

turbulent  reference  heat  transfer  rate, 
lbm/ft2-sec 

CFREFL 

Cf,  e,  ref.  I. 

laminar  reference  skin  friction  coefficient 

CFRKFT 

Cf,  o.  ref.  T 

turbulent  reference  skin  friction  coefficient 

A YE REF 

l\v,  ref 

reference  wall  enthalpy,  ft2/sec2 

20)  ROHMER 

RORMl'R  calculates  the  following  reference  properties: 
(a)  Viscosity  parameter 

1/14 


(b)  DensHy-viscotity  product 


Input 

POPS1.X 

P/PSL. 

local  pressure  along  streamline,  atm 

OMKGS 

stagnation  viscosity-temperature  ratio, 
slug-ft-sec-  ‘h 

OMEGE 

“e 

edge  viscosity-temperature  ratio,  s lug/ft- sec- °R 

OMEGW 

ww 

wall  viscosity-temperature  ratio,  slug/ft-sec- *R 

Output 

OMEGR 

r 

reference  viscosity-temperature  ratio, 
slug/ft-sec-  #R 

ROMUR 

P  P 
rrrr 

reference  density-viscosity  product,  slug'Vft^-sec 

ZTR 

ZTr 

reference  compressibility-temperature 
product,  ®R 

AYER 

V 

0  9 

reference  enthalpy,  ft'Vsec 

S1GR 

CTr 

reference  Prandtl  number 

21)  SiRCH  (Y,  X,  YA,  XA,  N,  IO) 

SIRCH  uses  TAB  to  perform  single  interpolation  on  a  table  of  y  =  f  (x). 
Input 


X 

X 

independent  variable 

YA 

dependent  array 

XA 

independent  array 

N 

number  of  values  in  XA  and  YA  arrays 

IO 

order  of  interpolation 

Output 

Y 

y 

dependent  variable 

r>8 


22)  SOM  EGA  {ENTH,  °OPSLX,  ZT,  OMEGA) 

SOMEGA  calculates  the  viscosity-temperature  ratio  w  and  the  compressibility- 
temperature  product  ZT  as  functions  of  enthalpy  and  pressure. 

(a)  Convert  enthalpy  to  Btu/lbm 

i(Btu/lbm)  =  25301  rift2 /sec  ) 

(b)  Calculate  ZT 
For  i  £  180 

ZT  =  i/.  432 

For  i  >  180 

ZT  =  f(log1Q  P/PSL.  i)  [TABLE2] 


Calculate  w 

ZT  <111  o>  =  14.454  x  10 


111  <  ZT  <  2000  u 


1/2 

/1.8  ZT  \  /  14.454  jr 

'  200  /  \1.8  ZT+  200  /  14‘ 


x  10 


-10 


ZT  >  2000  u  =  f(log10  P/PgL,  ZT)  (TABL18) 


(c)  Convert  ZT  and  u  to  *R 

ZTtR  =  1.8  ZT.k 


^slug  _  u 
ft-sec-*R  1.8 

Input 

ENTH  i  enthalpy 

POPSLX  P/PsL  local  pressure  along  streamline,  atm 


Output 


ZT 


compressibility-temperature  product,  *R 
viscosity-temperature  ratio,  slug/ft-sec-  *R 


ZT 

OMKGA 

23)  SON K NT  (FNTH,  SONIC) 

SONENT  calculates  the  speed  of  sound  as  a  function  of  free  stream  enthalpy.  The 
result  is  obtained  by  linear  interpolation  on  tabular  arrays  constructed  from  the 
equation: 


where 

T„  =  f<icc/RTo> 
=  f<T00) 


Input 

ENTH 

‘oo 

2 

free  stream  enthalpy,  ft  /sec 

Output 

SONIC 

acc 

speed  of  sound,  ft/sec 

24)  STACON  (PSTPNF,  AYEST,  VELP,  ACH,  ACHN,  H,  PINF,  VEL) 

STACON  calculates  the  stagnation-free  stream  pressure  ratio,  the  stagnation 
enthalpy,  and  the  modified  Newtonian  pressure  gradient. 

(a)  Stagnation  enthalpy 


(b)  Pressure  ratio  Pq/P^ 


f.O 


3.5 


(c)  Modified  Newtonian  pressure  gradient 


Interpolate  to  find  the  empirical  correction  factor  to  the  modified  Newtonian 
velocity  gradientas  a  function  of  MwcosA  (Swept  cylinder)  or  (Hemisphere). 

UMN  =  <UMN>  *  (CORRECTION  FACTOR)5 

For  the  correction  factor  see  Figure  4. 


(il 


Figure  4:  CORRECTION  FACTOR  FOR  THE  VELOCITY  GRADIENT  AT  A 

HEMISPHERE  STAGNATION  POINT  OR  CYLINDER  STAGNATION 
LINE 


Ii2 


Input 


ACM 

free  stream  Mach  number 

ACHN 

MW 

normal  Mach  number 

H 

H 

total  enthalpy 

PINT 

P=c 

free  stream  pressure,  lb/ft^ 

VEL 

Uoc 

free  stream  velocity,  ft/sec 

Output 

PSTPNF 

VP- 

stagnation  to  free  stream  pressure  ratio 

AYEST 

iST 

stagnation  enthalpy,  ft2/sec2 

YELP 

UMX 

modified  Newtonian  velocity  gradient,  l/sec 

25)  STREAM  (Program  C) 

STREAM  calculates  the  streamline  coordinate  x  as  a  function  of  the  geometry 
coordinates  |  and  tj  .  A  complete  description  of  the  geometry  and  streamline  pattern 
are  provided  by  the  coordinates  |  and  rj,  0e(|,b),  bmax(|)  and  bmaxU). 

0e(|.  r l)  must  be  provided  for  the  grid  of  |  and  n  that  lies  within  the  geometry 
boundaries.  With  0max(|)  known  along  the  upper  rj  geometry  boundary  the  complete 
streamline  pattern  on  the  body  is  defined.  bmax(!)  specifies  the  upper  rj  boundary 
thus  defining  the  geometry. 

The  coordinates  (If,  bf)  determine  a  point  on  a  particular  streamline.  The 
streamline  calculation  begins  at  (If,  bf)  and  proceeds  upstream  until  the  geometry 
boundary  is  reached.  The  calculation  of  the  streamline  is  done  by  the  following 
procedure; 

Let  (!x)j  =  If,  (bx)  -  bf,  and  (x)j  =  0.  Linear  interpolation  is  performed  on 
0e(|,b)  and  0max(|)  if  necessary,  to  determine  (0x)j.  Then 

(d”x)1  =  dxsin(9x), 

(dL)j '  4xcos  (ex)t 


03 


Then  the  coordinates  of  the  new  point  are: 


K), 

N.2  *  K),  -  K), 
h)2 " 


Continuing  this  procedure,  the  general  equations  are. 


The  streamline  calculation  is  completed  by  rearranging  the  £x,  Vx  and  x  arrays  such 
that  fjx,  17  and  x  are  measured  from  the  upstream  end  of  the  streamline. 

The  subroutine  STREAM  assumes  zero  streamline  divergence  due  to  body 
geometry 
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Input 


M  numl>er  of  t,  coordinates 

N  number  of  tj  coordinates 


XS1F 

«f 

final  streamline  £  coordinate,  ft 

ETAF 

nf 

final  streamline  of  n  coordinate,  ft 

DX 

dx 

incremental  streamline  distance,  ft 

XSI 

4 

geometry  coordinate,  ft 

ETA 

TJ 

geometry  coordinate,  ft 

ETAMAX 

'max 

body  edge  coordinate,  ft 

TH8MX 

^max 

streamline  angle  at  body  edge,  degrees 

THETAE 

9e 

streamline  angle,  degrees 

Output 

II 

number  of  points  calculated  along  streamline 

XS1X 

*x 

streamline  coordinate,  ft 

ETAX 

^x 

streamline  coordinate,  ft 

X 

X 

streamline  coordinate,  ft 

Til  ETAX 

*x 

streamline  angle,  degrees 

XF 

xf 

final  x-coordinate  on  streamline,  ft 

26)  TABLE2  (AYES,  PEDGEL,  TS,  CHECK) 

TABLE2  obtains  by  table  interpolation  the  compressibility-temperature  product 
as  a  function  of  pressure  and  enthalpy. 


Input 

AYES 

i 

enthalpy,  Btu/lbm 

PPLOC 

loSlO  p/pSL 

log  of  ratio  of  pressure  to  sea  level  pressure, 
pressure  in  atm 

Output 


TS 

ZT 

compressibility-temperature  product,  *R 

CHECK 

program  check  to  determine  success  of  interpolation 
(if  the  interpolation  is  not  successful,  an  error 
statement  is  printed) 

27) 

TABLE6  (TVV,  PPLOG.  AYEW) 

TABLE6  interpolates  to  find  wall  enthalpy  as  a  function  of  pressure  and  wall 
temperature  tables. 

Input 

TW 

T,v 

wall  temperature,  *K 

PPLOG 

loglO  p/PSL 

log  of  ratio  of  pressure  to  sea  level  pressure, 
pressure  in  atm 

Output 

AYEW 

*w 

wall  enthalpy,  Btu/lbm 

28) 

TABLE7  (AYEE.  PPLOG, 

AYEAYE) 

TABLE7  interpolates  on  a  table  of  pressure  and  edge  enthalpy  to  find  the  dissocia¬ 
tion  on  reaction  enthalpy  ratio  i£)/ie. 

Input 

AYEE 

*e 

edge  enthalpy,  ft2/see2 

PPLOG 

l°g  10  p/pSL 

log  of  ratio  of  pressure  to  sea  level  pressure, 
pressure  in  atm 

Output 

AYEAYE 

iD^e 

dissociation  reaction  enthalpy  ratio 

29) 

TABL14 

(ZTW,  S1GT) 

TABL14  obtains  by  table  interpolation  the  partial  Prandtl  number  as  a  function  of 
the  compressibility-temperature  product. 


Z  T 


compressibility-temperature  product,  ’ll 


Input 
7.TW 
Output 

SIGT  o  Prandtl  number 

30)  TABUS 

TABLlS  determines  by  table  interpolation  the  viscosity-temperature  ratio  as  a 
function  of  the  pressure  and  compressibility-temperature  product. 


Input 

ZT 

ZT 

compressibility-temperature  product.  °K 

PPLOG 

>o«10  P/PSI. 

log  of  ratio  of  pressure  to  sea  level  pressure, 
pressure  in  atm 

Output 

OMEGA 

U1 

viscosity-temperature  ratio,  slug/ft-sec-*K 

CHECK 

program  check  to  determine  success  of  interpolation 

(if  interpolation  is  unsuccessful,  an  error  comment 
is  printed) 


31)  TABL19 

TABL19  performs  table  interpolation  to  determine  the  compressibility-temperature 
product  as  a  function  of  pressure  and  viscosity-temperature  ratio.  The  table  is  com¬ 
posed  of  data  in  TABLK2  and  TABJ.1S. 

Input 

OMEGA  uj  viscosity-temperature  ratio,  slug/ft-sec -  °K 

FPLOG  l°fiio  log  ra^°  °f  Pressure  to  sea  level  pressure, 

pressure  in  atm 

Output 

ZT  ZT 


compressibility-temperature  product,  "K 


32)  TABL20 


TABL20  performs  table  interpolation  to  find  i/RT0  as  a  function  of  pressure  and 
compressibility-temperature  product. 


Input 

ZT  ZT 

compressibility-temperature  product,  0  K 

PPLOG  log10  P/PSL 

log  of  ratio  of  pressure  to  sea  level  pressure, 
pressure  in  atm 

Output 

AYERTO  i/RT0 

dimensionless  enthalpy 

33)  TBLP  (XT,  YT,  X,  NTAB, 

N> 

TBLP  is  a  linear  single  interpolation  routine.  To  save  search  time  in  finding  the 
location  of  the  independent  x  value  in  the  x  table,  the  search  is  begun  at  the  previously 
loca.ed  x  value. 

Input 

XT 

independent  variable  table 

YT 

dependent  variable  table 

X 

independent  variable 

NTAB 

number  of  values  in  the  x  table 

N 

location  in  x  table  at  which  the  search  is  begun. 

This  value  should  be  set  =  1  upon  the  first  entry 
to  this  routine.  Thereafter,  it  stores  the  location 
of  the  previous  search. 
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0.  INPUT-OUTPUT  DESCRIPTION 


a.  Data  Input  Preparation 


Data  are  input  to  this  program  via  punched  cards.  The  purpose  of  this  section  is 
to  deline  the  required  input  and  the  form  it  takes  on  the  cards.  Input  sheets  for  each 
option  or  section  of  the  program  are  shown.  Use  of  this  type  of  form  greatly  reduces 
the  amount  of  effort  required  to  obtain  results  from  the  program.  Data  may  be  key 
punched  directly  from  the  forms. 


1)  AXISYMMETRIC  OR  TWO-DIMENSIONAL  BODIES  (Option  D-l) 


Card 

Columns 

Format 

Description 

1 

1-6 

A6 

Type  of  case:  FLIGHT.  WT1.  or  WT2. 

This  card  indicates  the  type  of  free-stream 
conditions  to  be  input. 

1 

7-10 

— 

Leave  blank. 

1 

11-72 

10A6 

Title  or  some  description  of  case. 

Card  -2  will  consist  of  only  one  of  the  following  types  of  input:  Flight.  Wind 
Tunnel  1  (WT1),  or  Wind  'Funnel  2  (WT2).  Wind  'runnel  1  case  should  be  used  for 
tunnels  operating  in  the  ideal  gas  region,  while  Wind  Tunnel  2  ease  should  be  used  for 
high  energy  tunnels  such  as  shock  tubes. 

FLIGHT 

2 

1-10 

F10 

Flight  altitude,  ft. 

2 

11-20 

no 

Flight  velocity,  ft/see. 

WT1 

2 

1-10 

F10 

Free-stream  Mach  number. 

2 

11-20 

F10 

Free-stream  static  temperature.  R. 

2 

21-00 

no 

Free-stream  static  pressure,  lb/ ft-. 

WT2 

2 

1-10 

no 

Free-stream  velocity,  ft/see. 

2 

11-20 

F10 

•>  •) 

Total  enthalpy,  ff  ,'see“. 

Card 

Columns 

Format  Description 

2 

21-30 

F10 

Ratio  of  stagnation  pressure  to  sea  level 
pressure. 

2 

41-50 

F10 

.  2 

Free-stream  dynamic  pressure,  lb/ft  . 

3 

1-10 

F10 

• 

Radius  of  60°  swept  infinite  cylinder  for 
the  turbulent  reference  heating  rate,  ft. 

3 

11-20 

F10 

Radius  of  hemisphere  for  the  laminar 
reference  heating  rate,  ft. 

3 

21-30 

F10 

Reference  wall  temperature,  °R. 

4 

1-14 

A6 

Contains  the  word  AXISYMMETRIC  or 
TWODIMENSIONAL. 

5 

1-5 

15 

M  is  the  number  of  values  in  the  geometry 
tables  y(£)  and  £  .  3  <  M  <l  50. 

6 

— — 

8F10 

This  table  contains  M  values  of  £  ,  ft. 

7 

- - 

8F10 

This  table  contains  M  values  of  y(£) 
corresponding  to  £  in  the  preceding  table, 
ft. 

Note: 

A  minimum  of  three  points  must  be  input  into 
the  preceding  two  tables.  If  the  body  is  very 
slender,  more  points  should  be  used. 

8 

1-10 

F10 

Initial  x  location,  point  at  which 
calculation  begins,  ft. 

8 

11-20 

F10 

Increment  in  x,  ft.  dx  >  (x^  -  xo)/750. 

8 

21-30 

F10 

Edge  velocity  at  Xj  ,  ft/sec. 

9 

1-5 

15 

ITC  is  the  number  of  values  in  the  wall 

temperature  versus  surface  distance 
tables.  2£ITC£50. 

If  ITC^3  then  the  wall  temperature 
must  be  nonisothermal . 
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Card 

Columns 

Format 

Description 

10 

— 

8F10 

This  table  contains  ITC  values  of  surface 
distance,  ft. 

11 

8F10 

This  table  contains  ITC  values  of  wall  or 

surface  temperature  corresponding  to  the 

preceding  distance  table,  R.  If 

T  =  constant  then  enter  two  values, 
w 

12 

1-5 

15 

KpQ  is  the  number  of  printout  locations 
desired.  Kp0  <  500. 

12 

11-20 

F10 

Transition  Reynolds  number,  see 
definition  of  Rr  q  in  Program  A,  labeled 
QTRANS  in  section  2.  b.  18. 

13 

— 

8F10 

This  table  contains  Kpo  locations  at 
which  printout  is  desired. 

2)  HEMISPHERE  CALCULATION  (Option  D-2) 


Card 

1 

Columns 

1 

Format 

Description 

1 

2 

O 

I 

i  Same  as  Option  D-l 

o 

4 

i 

1-10 

A6 

Contains  title  of  program  -  HEMISPHERE. 

5 

1-10 

F10 

Hemisphere  radius,  ft. 

6 

1-5 

15 

ITC  is  the  number  of  points  in  the  Tw 
versus  x  tables.  2  < ITC  <50. 

If  ITC  >  3  then  the  wall  temperature  must 
be  nonisothermal. 

7 

— 

8F10 

This  is  a  table  of  ITC  values  of  x 
where  Tw  are  input,  ft. 

8 

8F10 

This  is  a  table  of  ITC  values  of  Tw 
which  correspond  to  locations  given  in  the 
preceding  x  table,  or  if  Tw  =  constant, 
enter  two  values  of  Tw. 

9 

1-5 

15 

Kpo  is  the  number  of  printout  locations 
desired.  KpO<  500. 

9 

11-20 

F10 

Transition  Reynolds  number,  see 
definition  of  Rr,Q  in  Program  A, 
labeled  QTRANS  in  section  2.  b.  18. 

10 

— 

8F10 

This  table  contains  Kpo  surface  distance 

locations  around  the  hemisphere  at  which 
printout  is  desired. 
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TABLE  (ITC  VALUES) 


3)  SWEPT  INFINITE  CYLINDER  (Option  D-3) 


Card  Columns  Format 


Description 


Same  as  Option  D-l 


1-8 


A6 


This  card  contains  the  word  CYLINDER  to 
allow  the  program  to  select  the  D-3 
options. 


1-10 


F10 


Cylinder  radius,  ft. 


11-20 


F10  Cylinder  sweep  angle,  measured  trom  a 

line  normal  to  the  flow.  0  <  A  <  90  . 


1-10 


11-20 


F10  is  the  abscissa  at  which  the  streamline 

calculation  begins,  ft.  £<■  must  be  large 
enough  so  that  the  streamline  has  very 
nearly  zero  slope  at  the  stagnation  line. 
Cases  have  been  run  by  the  authors  using 
=  10  and  20  ft  for  A  -10 3  and  GO  . 

F10  rjj-  is  the  ordi  ate  at  which  the  streamline 

calculation  begins,  ft. 


21-30 


1-5 


F10 


15 


dx  is  the  increment  in  distance  along  the 
streamline  used  in  the  boundary  layer 
calculations  (Program  A),  ft. 


d*>  Uf 


i  •>  l  /'l 

/GOO 


ITC  is  the  number  of  values  in  Tv,  and  p 
tables.  2  <  ITC  1  30.  If  ITC  >  3  the  wall 
temperature  is  nonisothennal. 


8F10  This  table  contains  ITC  values  of  p  . 

where  rj  is  the  distance  around  the  cylinder, 
ft. 


8F10  This  table  contains  ITC  values  of  T 

corresponding  to  the  p  table, 3 R. 


w 


Card 


Columns 


Format 


Description 


10  1-5 

10  11-20 


10  21-30 


11 


15  Kpo  is  the  number  of  x  printouts  desired 

along  the  streamline.  K?q  <  500. 

F10  Transition  Reynolds  number,  see 

definition  of  Rr,Q  in  Program  A,  labeled 
QTRANS  in  section  2.  b.  18. 

F10  At*  is  an  increment  used  in  Program  C  to 

locate  the  starting  location  of  second 
streamline,  ft. 

If  At*  is  small,  in  some  cases  the 
resolution  between  streamlines  near  the 
streamlines  origin  or  apex  may  be  very 
poor.  The  poor  resolution  leads  to  heat 
transfer  distributions  that  are  not  smooth. 
In  the  event  this  occurs  At*  should  be 
increased,  ft. 

8F10  XpQ  is  a  table  containing  locations  along 

the  streamline  at  which  printout  is 
desired. 
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4)  SHARP  DELTA  WING  (Option  D-4) 


Cards  Columns  Format 

1  I 

2  >  Same  as  Option  D-l 

3  I 

4  1-6  A-6 

5  1-5  15 


6  1-10  FIG 

6  11-20  F10 

6  21-30  F10 

6  31-40  F10 

6  41-50  F10 


7  1-10  F10 


Description 


Contains  the  word  DELTA. 

ITC  is  an  indicator  for  wall  temperature 
gradient.  If  the  wall  is  at  constant 
temperature  input  a  1  and  one  value  of 
wall  temperature  in  table. 

Angle  of  attack,  degrees. 

Wing  sweep  angle  measured  from  a  line 
normal  to  the  flow,  degrees. 

is  the  final  or  end  coordinate  along  the 
streamline  at  which  heating  data  is 
desired,  ft. 

rj£  is  the  final  or  end  coordinate  along  the 
streamline  at  which  heating  data  is  desired, 
ft. 


xj  is  the  point  along  the  streamline  at  which 
the  calculations  begin  and  may  have  to  be 
greater  than  zero  as  indicated  in  Figure  5. 
The  resulting  change  in  streamline  curva¬ 
ture,  due  to  interpolation,  near  the  apex 
produces  a  discontinuity  in  the  heating 
calculations,  ft. 


is  the  increment  in  £  which  the  pro¬ 
gram  uses  to  setup  an  ( £ . 77)  array  of 
streamline  angles,  velocity,  and  pressures. 


ft. 


A?  > 


if 

50 


si 


Figure  5:  STREAMLINES  NEAR  THE  APEX  OF  A  SHARP  DELTA  WING 


Card 


Columns 


Format 


Description 


11-20 


F10  At)  is  the  increment  in  17  the  program 

uses  to  setup  an  (£,tj )  array  of  streamline 
angles,  velocity  and  pressures,  ft. 


7  21-30  F10  dx  is  the  increment  in  distance  along  the 

streamline  at  which  calculations  will  be 
performed,  ft.  dx  >  (x*-  -  Xj)/750. 

Note:  If  Txv  is  a  constant  value,  enter  one  value  into 
card  11.  and  delete  cards  S,  9,  and  10. 


8  1-5  15  M  is  the  number  of  points  in  the  £  array. 

2  <  M  <  50. 

8  6-10  15  N  is  the  number  of  points  in  the  n  array. 

2  <  N  <50. 

9  -  8F10  This  is  a  table  of  M  values  of  t  where 

wall  temperature  will  be  input,  if  ITC  >  1. 

10  -  8F10  This  is  a  table  of  N  values  of  n  where 

wall  temperature  will  be  input  if  ITC  >  1. 

11  -  8F10  This  is  a  table  containing  (N)  •  (M)  values 

of  wall  temperature,  input  so  i  varies  with 
each  jj  . 

12  1-5  15  Kp()  is  the  number  of  x  printouts  desired 

along  the  streamline.  Kp(,  <  500. 

12  11-20  F10  Transition  Remolds  number,  sic 

definition  of  Rr  q  in  Program  A.  labeled 
QTRANS  in  section  2.  b.  18. 

12  21-30  F10  See  card  10,  swept  infinite  cylinder. 

option  D-3. 


8F10  Xpq  is  a  table  of  Kpo  locations  along  the 

streamline  where  printout  is  desired. 


5)  STREAMLINE  CALCULATIONS  (PROGRAM  C) 


Card 

1 

Columns 

Format 

Description 

1 

2 

O 

i  Same  as  Option  D-l 

i) 

4 

i 

1-6 

A6 

Must  contain  the  word  STREAM. 

5 

1-5 

15 

M  is  the  number  of  points  in  the  array 

2  <  M  <  50. 

5 

6-10 

15 

N  is  the  number  of  points  in  the  array  n 

2  <  N  <50. 

5 

11-15 

15 

ITC  is  an  indicator  used  to  test  whether 
the  wall  temperature  is  a  constant  or  a. 
variable.  If  Tw  is  a  constant  ITC  =  1  and 
enter  one  value  of  T  into  card  16  and 
delete  cards  13.  14.  and  15. 

6 

1-10 

F10 

4  f  is  £  at  the  final  or  end  point,  ft. 

6 

11-20 

F10 

Tjr  is  tj  at  the  final  or  end  point,  ft. 

6 

21-30 

F10 

dx  is  the  increment  in  distance  along  the 
streamline,  ft. 

1.2  „2\1/2 

,  +  "f  > 

dx  4  750 

6 

31-40 

F10 

Xj  is  the  point  along  the  streamline  at 
which  the  calculations  will  begin,  ft. 

6 

41-50 

F10 

uc  Xj  is  the  edge  velocity  at  the  initial  or 
start  of  calculations,  ft/sec. 

7 

— 

8F10 

This  is  a  table  of  M  values  of  £ 

2  <  M  <  50.  8q.  P/P0  and  Tw  will  be 

input  as  functions  of  £  and  n  .  ft. 


Card 


Columns 


Format 


Description 


8  8F10 

9  8F10 

10  8F10 

11  8F10 

12  8F10 

13  1-5  15 

13  6-10  15 

14  8F10 

15  8F10 

16  8F10 

17  1-5  15 


This  is  a  table  of  N  values  of  tj 
2  <  M  <  50.  6e,  P/PQ  and  Tw  will  be 

input  as  functions  of  £  and  rj ,  ft. 

This  is  a  table  of  M  values  of  nmax  as  a 
function  of  £ .  TJmax  defines  the  upper  or 
outer  boundary  of  the  body,  ft. 

This  is  a  table  of  M  values  of  the  stream¬ 
line  angle  0max  along  the  boundary  as  a 
function  of  £ ,  deg. 

This  is  a  table  of  (M)  .  (N)  values  of  the 
local  streamline  angle  8e  and  should  be 
input  along  lines  of  constant  tj  and  £  , 
deg. 

This  is  a  table  of  (M)  •  (N)  values  the 
local  pressure  divided  by  the  stagnation 
pressure  and  should  be  input  along  lines  of 
constant  tj  and  varying  £ . 

MT  is  the  number  of  £  values  entered  into 
the  table  Tw  =  f(£,r/)-  Also  see  ITC  above. 

NT  is  the  number  of  tj  values  entered  into 
the  table  Tw  =  f(  £  ,  tj  ).  Also  see  ITC  above. 

This  is  a  table  of  MT  values  of  £  at  which 
Tw  will  be  input.  Also  see  ITC  above. 

This  is  a  table  of  NT  values  of  tj  at  which 
Tw  will  be  input.  Also  see  ITC  above. 

This  is  a  table  of  (MT)  •  (NT)  values  of 
Tw  as  a  function  £  and  tj  ,  °  R.  Tw 
should  be  input  as  a  function  of  £  at 
constant  tj  . 

Kpo  is  the  number  of  x  locations  along 
the  streamline  at  which  printout  is  desired. 
KpQ  <  500. 
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Card 

Columns 

Format 

Description 

17 

11-20 

FIG 

Transition  Reynolds  number,  see  defini¬ 
tion  of  Rr.Q  iii  Program  A.  labeled 
QTRANS  in  section  2.  b.  18. 

17 

2i-:;o 

F10 

See  card  10.  Swept  infinite  cylinder, 
option  D-2. 

18 

1-10 

FI  0 

xpo  is  a  table  of  Kpo  locations  along  the 
streamline  at  which  printout  is  desired. 

table  <n  values 


6)  REFERENCE  CALCULATIONS  (PROGRAM  REF) 


Card  Columns  Format  Descriptor. 

1  \ 

2  )  Same  as  Option  D-l 

3  ] 

4  1-10  A6  Must  contain  the  word  REFERENCE . 


7)  BOUNDARY  LAYER  CALCULATIONS  (PROGRAM  B) 


Card  Columns  Format  Description 

1  | 

2  >  Same  as  Option  D-l 

3  ) 

4  1-6  A6  Must  contain  the  word  FLOW. 

5  1-5  15  IS  is  the  number  of  values  in  the  P/P0, 

8e  ,  A/Aj  and  t/t[  tables  as  a  function  of 
streamline  distance.  2  <  IS  <  50. 

6  8F10  This  is  a  table  containing  IS  number  of  x 

locations,  ft,  along  the  streamline  at  which 
P/P0,  0e,  A/Aj  and  r/rj  will  be  input. 

7  8F10  This  is  a  table  containing  IS  values  of 

local  pressure  each  normalized  by  the 
model  or  vehicle  stagnation  pressure.  The 
P/PQ  correspond  to  the  above  x  table. 

8  8F10  This  table  contains  IS  values  of  the  local 

streamline  angle  relative  to  the  axis  of 
symmetry,  degrees.  The  values  of  8Q 
correspond  to  the  above  x  table,  degrees. 

9  8F10  This  table  contains  IS  values  of  the  stream¬ 

line  divergence  parameter  A/Aj .  The 
values  of  A/A^  correspond  to  the  above  x 
table. 

10  —  8F10  This  table  contains  IS  values  of  the  body 

shock  geometry  parameter  r/r^.  The 
values  of  r/ r j  correspond  to  the  above  x 
tabic. 

11  1-10  F10  xj  is  point  at  which  the  calculation  will 

begin,  ft. 


Card 


Column 


Format 


Description 


11 

11-20 

F10 

d>:  »s  the  increment  in  distance  along  the 
streamline,  ft 

,  N  FINAL  XI 

dx  > 

/o0 

11 

21-30 

F10 

ue  N  is  the  edge  velocity  a  ft  sec. 

12 

1-5 

15 

ITC  is  the  number  of  values  in  the  T 
versus  x  table.  2  <  ITC  <  50. 

If  Tu  =  constant  then  ITC  -  2. 

16 

— 

8F10 

This  is  a  table  of  x locations  (ITC)  along 
the  streamline  where  Tw  will  be  input. 

14 

8F10 

This  is  a  table  of  ITC  values  of  Tu.  along 
the  streamline  corresponding  to  the  x 
table.  If  ITC  =  2  enter  only  two  values  of 
TW,  *R. 

la 

1-5 

15 

Kpo  is  the  number  of  x  locations  along  the 
streamline  at  which  printout  is  desired. 

KpQ  <  500. 

15 

11-20 

F10 

Transition  Reynolds  number,  see 
definition  of  Rr  q  in  Program  A.  labeled 
QTRAXS  in  section  2.  b.  18. 

1G 

— 

8  F10 

This  is  a  table  of  Kpf)  locations  along  the 

streamline  where  printout  is  desired. 


96 


TABLE  (!S  VALUES) 


A/A.  TABLE  (IS  VALUES) 


8)  HEAT  TRANSFER  CALCULATIONS  (PROGRAM  A) 


Card 

1  \ 

Columns 

Format 

Description 

1  | 

2  \ 

o  1 

Same  as  Option  D-l 

<>  1 

4 

1-6 

A6 

Must  contain  the  word  QTRANS. 

o 

1-5 

15 

IS  is  the  number  of  values  in  each  of  the 
following  tables:  x.  P/PsL’  ue- 
A/Aj  and  r/rj .  2  <  IS  <  50. 

6 

1-10 

F10 

Xj  is  the  initial  x  location  aiong  the  stream 
line  at  which  the  calculations  will  start,  ft. 

6 

11-20 

no 

dx  is  the  increment  in  distance  along  the 
streamline,  (ft. ) 

,  XFINAL-XI 

dx  > - — — - 

i  oO 

6 

21-30 

F10 

ue..\j  ‘s  boundary  layer  edge  velocity 

atXj,  ft/sec. 

7 

" 

8F10 

This  is  a  table  containing  x  locations  (IS) 
along  the  streamline  at  which  P/PsL’  ue  • 
9e  ,  A/Aj  and  r/rj  will  be  specified,  ft. 

8 

8F10 

This  is  a  table  containing  IS  values  of  the 
local  pressure  normalized  with  the  sea 
level  pressure,  and  at  locations  correspond 
ing  to  the  x  table. 

9 

6F10 

This  is  a  table  containing  IS  values  of  the 
boundary  layer  edge  velocity  correspond¬ 
ing  to  locations  specified  in  the  x  table, 
ft/ sec. 

10 

— 

8F10 

This  is  a  table  containing  IS  values  of  the 

local  streamline  angle  f?c  relative  to  the 
axis  of  symmetry.  The  angles  correspond 
to  x  locations  specified  in  the  x  table. 


deg. 
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Card 

Columns 

Format 

Description 

11 

8F10 

This  is  a  table  containing  IS  values  of  the 
streamline  divergence  parameter  A/A.- 
and  correspond  to  x  locations  specified  in 
the  x  table. 

12 

8F10 

This  is  a  table  containing  IS  values  of  the 
body-shock  parameter  r/r^  and  correspond 
to  x  locations  specified  in  the  x  table. 

13 

1-5 

15 

ITC  is  the  number  of  values  in  the  T„. 
versus  x  table. 

2  <  ITC  <  50 

If  Tw  is  a  constant,  ITC  -  2.  Input  two 
values  of  Tw  into  card  15. 

14 

— 

8F10 

This  is  a  table  containing  ITC  values  of 
x  along  the  streamline,  ft. 

15 

- - - 

8F10 

This  is  a  table  of  ITC  values  of  Tw  along 
the  streamline  and  corresponding  to  the 
x  locations  in  the  preceding  table.  (c  R) 

16 

1-5 

15 

KpQ  is  the  number  of  x  printout  locations 
desired  along  the  streamline.  Kpt)  <  500. 

16 

11-20 

F10 

Transition  Reynolds  number,  see  defini¬ 
tion  of  Rr,Q  in  Program  A,  labeled 
QTRANS  in  section  2.  b.  18. 

16 

21-30 

F10 

BEQL  is  the  ratio  of  the  laminar  equi valent 
distance  parameter  to  the  distance  along 
the  streamline.  If  xj  >  0  then  a  value  of 
BEQL  is  required.  If  xj  =  0  input  BEQL=0. 

16 

31-40 

F10 

BEQT  is  the  turbulent  equivalent  distance 
parameter.  See  BEQL  above. 

17 

— 

8F10 

This  is  a  table  of  KPq  locations  along  the 

streamline  where  printout  is  desired. 
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TY 


b.  Output  Description 


free  stream  velocity,  ft/sec 

free  stream  Mach  number 

free  stream  enthalpy,  ft2/sec2 

free  stream  static  temperature,  °R 

2 

free  stream  static  pressure,  lb/ ft 
free  stream  dynamic  pressure,  lb/ft2 
(model  total  pressure)/ (sea  level  static  pressure) 

(model  total  pressure)/(free  stream  static  pressure) 

REFERENCE  CONDITIONS 

HO  %/(ia\v~i\v)  =  Vcp  "here  ho  is  the  hemisphere  stagnation 

point  heat  transfer  coefficient,  lbm/ft^-sec 

.  24HO  .24  ho/Cp,  Btu/ft2-sec-°  F. 

2 

QDOTO  hemisphere  stagnation  point  heating  rate,  Btu/ft  -sec 

HREF,L  ^REF,  l/(‘aw ,  L  "  Sv,  REF^  =  hREF,  L^CP  where  hREF.  L 

is  the  laminar  heat  transfer  coefficient  based  on  enthalpy  at 

the  stagnation  line  of  a  60°  swept  infinite  cylinder,  lbm/ 
f^-sec 

RHEM.REF  hemisphere  radius,  ft 

RCYL,  REF  cylinder  radius,  ft 

SWEEP  sweep  angle  of  the  swept  cylinder  used  for  the  turbulent 

reference  condition,  degrees 

.  24HREF,  L  ,24  hp^jrj?  j^/Cp 

QDOT,  REF,  I.  laminar  heating  rate  at  the  stagnation  line  of  a  60°  swept 

infinite  cylinder,  Btu/ft^-sec 


INITIAL  CONDITIONS 

VELOCITY 

MACH 

ENTHALPY 

TINF 

PINF 

QINF 

PO/PSL 

PO/PINF 
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HREF,T 


-  24HREF.  T 
QDOT,  REF,T 


AREF.T^aw.T  "  *w.REF)  ~  hREF,T/cP  where  hREFX  is 
the  turbulent  heat  transfer  coefficient  at  the  stagnation  line 
of  a  60°  swept  infinite  cylinder  based  on  enthalpy,  lbm/ft^-sec 

24  hREF,  t/  CF 

turbulent  heating  rate  at  the  stagnation  line  of  a  GO  swept 
infinite  cylinder,  Btu/ft“-sec-c  R 


IAW,  REF,  L 


laminar  adiabatic  wall  enthalpy  on  a  60°  swept  infinite 
cylinder,  ft^/sec2 


IAW,  REF,  T 

IW,  REF 


turbulent  adiabatic  wall  enthalpy  on  a  60°  swept  infinite 
cylinder,  ft2/sec2 

reference  wall  enthalpy,  ft^/sec^ 


TAUREF,  L 

CFINF,  REF,  L 

IW.O 


laminar  shear  stress  on  a  60;  swept  infinite  cylinder, 
lbj/ft2 

^Tref  i/(p«u«°)’  skin  friction  coefficient  based  on 

free  stream  conditions 

stagnation  wall  enthalpy,  ft^/sec2 


TAUREF,  T 

CFINF,  REF,  T 


turbulent  shear  stress  on  a  60°  swept  infinite  cylinder. 

Ibf/ft2 

^Tref  turbulent  skin  friction  coefficient  based  on 

free  stream  conditions 


TW,  REF 


wall  temperature  for  60c  swept  infinite  cylinder  heating 
calculations,  °  R 


GEOMETRIC  PARAMETERS 

SWEEP  CYL  cylinder  sweep  angle,  degrees 

SWEEP  DELTA  WING  delta  w  ing  sweep  angle,  degrees 
ALPHA  angle  of  attack,  degrees 

RCYL  cylinder  radius,  ft 

RHEMI  hemisphere  radius,  ft 


10.') 


STREAM  LINE 


XSI  coordinate,  ft 

ETA  coordinate,  ft 

THETAE  streamline  angle,  degrees 


BOUNDARY  LAYER  CALCULATIONS 


X 

ETA 

XSI 

BEQXL 

BEQX~ 

POPSL 

A  YE  AWL 

A  YE  AWT 

UE 

HL 

HT 

DODI 

N 

QLISO 

QTISO 

RORI 


streamline  distance,  ft 

ordinate  of  coordinate  system,  ft 

abscissa  of  coordinate  system,  ft 

(beq/x)L-  laminar  equivalent  distance  parameter 

(beq/x)j,  turbulent  equivalent  distance  parameter 

local  pressure  in  atmospheres 

2  2 

iaw>  l,  adiabatic  wall  enthalpy  in  laminar  flow,  ft  /sec 
iaw  X’  ac^abatic  wall  enthalpy  in  turbulent  flow,  ft  /sec" 
edge  velocity,  ft/sec 

laminar  heat  transfer  coefficient  divided  by  the  reference 
hemisphere  stagnation  point  heat  transfer  coefficient 

turbulent  heat  transfer  coefficient  divided  by  the  reference 
turbulent  heat  transfer  coefficient  at  the  stagnation  line  of  a 
60a  swept  infinite  cylinder 

A/A^,  total  streamline  divergence 

x/Za/Aj)  ^(A/Ajj/dx),  streamline  divergence  parameter 
laminar  isothermal  heat  transfer  rate,  Btu/ft^-sec-“  R 
turbulent  isothermal  heat  transfer  rate,  Btu/ft^-see-°  R 
r/rj,  streamline  divergence  due  to  body-shock  geometry 


10(5 


QDOTL 

QDOTT 

TW 

CFEL 

CFET 

ZTEX 

RRQ 

RRS 

OMEGAE 
AIEE 
A  YEW 
THETA 
MUO 

ZTR 

EBARL 

EBART 

JL 

XEQL 

FXQ 

FXS 

OMFGAR 


o 

laminar  heating  rate,  Btu/ff^-see-  R 

2 

turbulent  heating  rate,  Btu/ft“-see-"  R 
wall  temperature.  R 

laminar  skin  frietion  coefficient  \  based  on  dynamic 

>  pressure  at  boundary 

turbulent  skin  frietion  coefficient  )  layer  edge 

ZT  ,  compressibility  factor  times  the  edge  temperature,  UR 

heat  transfer  reference  Reynolds  number 

skin  frietion  reference  Reynolds  number 

2 

edge  viscosity  divided  by  Z'I'C.  lbj.-sec/ft“-  R 
2  2 

i  edge  enthalnv.  ft-/ see" 
e 

0  2 

1  ,  wall  enthalpy.  ft"/see“ 

local  streamline  angle,  degrees 

H  ,  absolute  viscosity  at  the  stagnation  reference  condition, 
lbj-scc/ft“ 

compressibility  factor-temperature  product  evaluated  at  the 
reference  enthalpy,  '  R 

laminar  crossflow  pressure  gradient  parameter 

turbulent  crossflow  pressure  gradient  parameter 

laminar  strcamwisc  pressure  gradient  profile  parameter 

laminar  heat  transfer  equivalent  distance  divided  by  the  local 
distance  along  the  streamline 

see  definition  of  F  __  in  QTRANS 

x.  Q 

see  definition  of  F  in  QTRANS 

*\«  u 

reference  viscosity  divided  by  reference  temperature. 
lbf-sec/ft2-  R 
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THETA L 


laminar  momentum  thickness,  ft 


THETAT  turbulent  momentum  thickness,  ft 

c.  Error  Statements 

ERROR  IN  CONTROL  CARD  (NTYPE  -  name) 

ERROR  IN  CONTR  C.  CARD  (NPROG  -  name) 

The  above  errors  are  caused  by  an  input  error  of  the  code  words  NTYPE  or 
NPROG.  "Name"  is  the  control  word  that  was  input.  This  error  will  terminate  the 
run. 


ALTITUDE  INPUT  TO  ATMOS  EXCEEDED  RANGE  -  ASSIGN  VALUE  OF 
2,300.000  FT. 

ALTITUDE  INPUT  TO  ATMOS  BELOW  ROUTINE  RANGE  -  ASSIGN  VALUE  OF 
0  FT. 

Input  altitude  out  of  range  of  ATMOS  subroutine.  FLIGHT  case  only. 

STREAMLINE  COORDINATE  DIMENSION  (750)  EXCEEDED  IN  SUBROUTINE 
AXED 

The  input  value  of  dx  is  too  small  for  an  axisymmetrie  or  two-dimensional  case. 
Also  printed  with  this  comment  is  the  x-value  at  which  the  array  dimension  was 
exceeded,  the  maximum  x  defined  by  the  geometry  input,  and  the  input  dx.  The  run 
is  terminated  by  this  error. 

DBTP  ERROR 

The  independent  variable  is  outside  the  range  of  the  program  tables.  The 
independent  variable  and  the  table  limits  are  printed.  The  program  will  use  the  table 
limit  and  proceed. 

STREAMLINE  CROSS  AT  XSI  = 

When  the  twc  divergence  calculation  streamlines  cross,  the  above  comment  is 
printed  with  the  value  of  the  4  coordinate  at  that  point.  The  program  moves  down¬ 
stream,  by  dx  increments,  until  the  streamlines  are  not  crossed  and  begins 
calculations  at  that  point. 


10.H 


CONVERGENCE  NOT  ESTABLISHED  IN  VELOCITY  CALCULATION 


The  velocity  iteration  routine  in  subroutine  FLOW  failed  to  converge  on  a 
solution  in  10  iterations.  The  program  proceeds  using  the  final  calculated  velocity. 

S1RCH  ERROR 

Table  interpolation  error  in  the  general  interpolation  routine  SIRCH. 

ERROR  IN  TABLE2  LOOK-UP 
ERROR  IN  TABU?  LOOK-UP 

An  error  occurred  during  table  interpolation  in  subroutine  TABLE2  or  TABL18. 

FREE  STREAM  ENTHALPY  EXCEEDS  TABLE  VALUES  IN  SON  ENT 

The  enthalpy  value  used  for  table  interpolation  in  subroutine  SON  ENT  is  outside 
the  table  limits. 

ERROR  IN  STREAMLINE  CALCULATION 

This  comment  is  printed  when  a  calculated  streamline  point  lies  outside  the 
.Tj)  grid  in  subroutine  STREAM. 
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4.  PROGRAMMING  INFORMATION 
a.  Program  listing 
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READ (5 *5010) <TWS( I >  *  I -1*1 TC)  24190120 
READ (5  *5050)KF»RTRANS*DETAS  24190121 
READ (5 *5010) (XP0( I ) *  I *1 *KF )  2419012  2 
CALL  CYLIND  24190123 
CALL  STREAM  24190124 
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b.  PrPr  Deck  Setup 


The  following  diagrams  indicate  the  deck  setup  and  the  core  storage  requirements 
for  the  Prpr  program.  Each  link  is  composed  of  the  subroutine  decks  (symbolic  or 
object)  listed  in  the  allocation  diagram.  All  control  cards  are  indicated  by  a  $  in 
column  1  and  must  be  punched  exactly  as  shown  in  the  deck  setup  diagram.  All  items 
in  parenthesis  indicate  subroutine  or  data  card  decks. 
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SYSTEM 
Link  0 


__  00000 
LP1 _ 


LP2 

LP3 


Notes: 

1.  Storage  locations  are  on  actual  numbers. 

2.  Load  points  (LP)  are  indicated  as  follows: 

LP1  -  04076 
LP2  -  60005 
LP3  -  65623 

SCHEMATIC  DIAGRAM  OF  COMPUTER  STORAGE  ALLOCATION 
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DECK  SETUP  INSTRUCTIONS  AS  2419  IBM  7094 
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The  "ontmits  of  each  link  are: 


lank  0 

AS2419  (main  program) 

TAB 

DTAB 

DBTP 

TBLP 

XTAB 

LOC 

XLOC 

Link  1 

DIVERG 

FIND 

STREAM 

Link  2 

W’ALLT  1 
W'ALLT  2 
BLKDTA 
DERV 
EBAR 
IW’ALL 
JAYELL 
MUZERO 
RORMUR 
SOM  EGA 
TABLE  2 
TABLE  6 
TABLE  7 
TABLE  14 
TABLE  18 
TABLE  19 
TABLE  20 


Link  3 

DELTA 

Link  4 

CY  LIN'D 
FLOW 

Link  5 

AXI2D 

HEMI 

Link  6 

QTRAN 

BEQI 

Link  7 

INTI  A  L 

REF 

SONENT 

STACON 

ATMOS 


SECTION  HI 


THE  TURBULENT  NONSIMILAR  BOUNDARY  LAYER  PROGRAM  (NSBL) 


1.  NOMENCLATURE 

Program 

Symbol 

Math 

Symbol 

Description 

Units 

CP 

CP 

specific  heat  of  air  at  constant 
pressure,  6006 

2  2 
ft  /sec  -°R 

6 

boundary  layer  thickness 

ft 

*  DELETA 

6(tj  ip ) 

increment  in  boundary  layer  parameter 

ft3/slug 

DELSTR 

6* 

boundary  layer  displacement  thickness 

ft 

*  DIIDX 

ai/ax 

derivative  of  total  enthalpy 

ft/ sec2 

♦  DHDY 

ai/ay 

derivative  of  total  enthalpy 

ft/sec2 

'  DMUDY 

dji/  dV 

derivative  of  viscosity 

lb-sec/ft3 

DQ 

dq/dy 

derivative  of  the  heating  rate 

3 

Btu/ft  -sec 

DRDX 

hi/  ax 

derivative  of  three-dimensional  flow 
parameter 

— 

DS 

AS 

increment  for  stagnation  region  problems 

ft 

*  DUDX 

3u/9x 

derivative  of  velocity 

l/sec 

*  DUDY 

9u/ 9y 

derivative  of  velocity 

1/sec 

*  DUDYL 

(9u/9y)e 

velocity  derivative  at  boundary  layer 
edge 

l/sec 

*  DVDY 

9v/9y 

derivative  of  velocity  normal  to  surface 

l/sec 

DX 

Ax 

increment  in  x-direction 

ft 

*  Those  symbols  marked  with  an  asterisk  are  arrays  (tables)  in  the  program. 
Throughout  this  document,  a  subscript  i  on  either  the  program  symbol  or  the 
equivalent  math  symbol  indicates  that  reference  is  being  made  to  the  i*h  position 
within  the  array. 


Program 

symbol 

Math 

symbol 

Description 

Units 

DY 

Ay 

increment  in  y-direction 

ft 

*  D2HDY 

a2i/ay2 

second  derivative  of  total  enthalpy 

w  2 

1/sec 

*  D2UDY 

a2u  /ay2 

second  derivative  of  velocity 

1  /sec-ft 

ENDX 

xf 

final  value  of  x 

ft 

*  ENTH 

i 

static  enthalpy 

f2,  2 
ft  /sec 

EPS 

€ 

eddy  viscosity 

lb-sec/ft2 

EPSY 

de/dy 

derivative  of  eddy  viscosity 

lb-sec/ft2 

*  ETA 

TJ 

similarity  parameter 

— 

*  H 

I 

total  enthalpy 

ft2 /sec2 

*  PR 

a 

Prandtl  number 

— 

*  PRESSURE 

P 

pressure 

lb/ft2 

*  QP 

q 

heating  rate 

,  2 

Btu/ft  -sec 

QW 

heating  rate  at  the  wall 

Btu/ft2-sec 

QQW 

QQW 

heating  rate  at  the  wall 

2 

Btu/ft  -sec 

R 

r 

three-dimensional  flow  parameter 

ft 

R 

universal  gas  constant,  1716 

ftVsec2-0  R 

*  RHO 

P 

density 

2  4 

lb-sec  /ft 

SI 

si 

initial  value  of  S  for  the  calculation  of  17 

ft 

STARTX 

X 

0 

initial  value  of  x 

ft 

*  t 

T 

temperature 

0  R 

*  TAU 

T 

shear  force 

lb/ft2 

TAULP 

tl 

local  laminar  shear  force 

lb/ft2 
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Program 

symbol 

Math 

symbol 

TAULW 

T  w,  L 

TAULY 

3TL/3y 

TAUTP 

rT 

TAUTW 

tk,T 

TAUTY 

3rT/3y 

TAUW 

Tw 

THETA 

0 

TURBS 

3t  /dy 

Va>' 

*  u 

u 

*  V 

V 

X 

X 

XI 

*1 

*  XMU 

V 

*  y 

3' 

Subscript  notation: 

Subscript 

e 

i 

I 

L 

o 

T 

w 


Description  Units 

/  2 

laminar  shear  force  at  the  wall  lb/ft 

3 

derivative  of  the  laminar  shear  force  lb/ft 

2 

local  turbulent  shear  force  Ib/ft 

/  2 

turbulent  shear  force  at  the  wall  lb/ft 

3 

derivative  of  the  turbulent  shear  force  lb/ft 

2 

shear  force  at  the  wall  lb/ft 

boundary  layer  momentum  thickness  ft 

ratio  of  the  turbulent  to  the  laminar 
shear  force  derivative 

velocity  in  the  x-direction  ft/sec 


velocity,  v’  in  Volume  I  (Appendix  C)  ft/sec 
coordinate  tangent  to  the  surface  ft 

initial  value  of  x  for  the  calculation  of  p  ft 

2 

absolute  viscosity  of  air  lb-sec/ft 

coordinate  normal  to  the  surface  ft 


evaluated  at  the  boundary  layer  edge 

evaluated  at  the  i^  position  within  an  array 

initial  value 

laminar 

initial  value 

turbulent 

evaluated  at  the  surface 
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2.  METHOD  OF  SOLUTION 
a.  Forward  Integration 


In  the  following  discussion,  the  equation  numbers  in  parenthesis  refer  to  an 
appropriate  equation  in  Appendix  C  of  Volume  I.  The  equations  solved  by  this  com¬ 
puter  program  are  the  equation  of  state  (C-3), 

P  -  P/(RT) 


Continuity  (C-17), 


V  —  - 

1 

9u . 

2u.l 

i 

i 

3v 

Ay 

_1  9P 
P  9x  + 


_a_ 

ay 


x-momentum  (C-6), 

ap  a  /  au  \1  vau 
ax  +  ay  ^  ay+  Tj j  u  ay 

and  Energy  (C-7): 


3u  _  l_ 
dx  pu 


ai_ 

dx 


J_  JL  ii  £L  *  • 

pu  3y  a  3y  +  4  U 


v  ai_ 
u  a> 


In  addition,  the  relationship  for  turbulent  shear  stress,  Equation  (C-39),  is 
"equired  for  the  computation. 


<V 


I  11 

a  9y 


Since  v  is  expressed  as  a  function  of  input  data  it  can  be  determined  explicitly  at 
each  point  in  the  boundary  layer  at  the  initial  or  start  position. 

The  y  derivatives  in  the  above  expressions  are  calculated  using  three-point 
central  differences  at  all  positions  except  at  y  =  Ay.  Since  the  turbulent  velocity 
and  enthalpy  profiles  change  very  rapidly  near  the  wall,  an  off-center  method  of 
calculation  of  derivatives  is  used  at  y  =  Ay. 


With  v  defined  and  with  u,  I  and  y  derivatives,  3u/3x  and  31/3 x  can 
now  be  determined.  With  the  3u/3x  and  31/ 3x  derivatives  determined  the  u  and  I 
profiles  at  the  next  station  (xQ^  +  Ax)  can  be  obtained  by  forward  intt  ation  using: 


The  y  derivatives,  v,  3u/3x,  and  3I/3x  calculations  move  out  from  the  wall  along  a 
line  normal  to  the  wall  until  a  specified  limit  is  reached.  The  profiles  are  then 
stepped  forward  and  the  calculations  repeated. 

At  each  point  in  the  boundary  layer,  a  similarity  parameter  i)  is  calculated  as 
indicated  below  in  subsection  b. ,  item  (4).  A  value  of  hmax  is  an  item  of  input  used 
to  limit  the  calculation  in  the  y-directions,  which  are  to  be  printed  out. 

Also  calculated  at  each  station  x  are  the  displacement  and  momentum  thickness, 
heating  rate,  and  the  laminar  and  turbulent  shear  stress  at  the  wall  using: 


1  The  notation  x0  will  be  used  throughout  this  section  to  denote  the  initial  value  of  x. 
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Sv 


1 

778 


dy 


r 


L 


djj_  9u 

9y  9y 


dy 


The  heating  rate  and  shear  at  the  wall  are  calculated  by  an  integration  over  the 
boundary  layer  rather  than  from  the  definitions  because  of  the  greater  accuracy 
which  can  be  obtained  with  the  equation. 

b.  Numerical  Method  of  Solution 


The  purpose  of  this  section  is  to  give  the  procedure  used  in  solving  the  partial 
differential  equations  defined  previously  including  the  required  numerical 
approximations. 

The  following  variables  and  tables  are  input.  Tables  are  indicated  as  functions 
of  some  independent  variable,  e.g.  P  =  f(x)  indicates  the  table  of  pressure  as  a 
function  of  x. 


u  =  f(y)  at  xQ 

I  «  f(y)  at  xq 

9u 

ay 

=  f(x) 

e 

P  -  f(x) 

\ 

»-* 

li 

£ 

r  -  f(x) 

i  =  f(x ' 
w 

In  addition,  the  following  quantities  must  be  specified: 
Ax  Ay  n. 


max 


f 


(1)  Calculation  of  tabular  values 


x.  -  x.  ,  x.  ,-x.l  r  , 

(P  P.)J_J1+(P..P.  )_*L_I - - - 

1+1  lxi+l"xi  1  1-1  xi-xi-lilxi+l"xi-l 


i  =  2,  3,  .  .  . ,  (NP-1) 


ap 

ax 


i=l 


=  0 


(2)  Values  obtained  from  tables  (linear  interpolation) 


P  =  f(x) 

(ij)  =fw 

y  e 

V!  =  f' 

r  -  f(x) 

‘w  =  f<X> 

“  f(x) 

Computation  at  initial  x  =  x^ 

y.  =  (i  -  1)  Ay 

i  =  1,  2,  3 . 200 

Wi  +  (l?)  *-v 

J  e 

i  =  (LIST  +  1),  ....  200 

\  !list 


LIST  =  number  of  input  values  in  u  and  I  tables. 
T.  =(I-u2/2)/Cp 

o  9  /9 

p.  =  (2.272  x  10  T.'  )/(T.  +  198.6) 

P.  =  P./(1716  T.) 
i  i  '  r 

Cp  -  6006  ft2/sec-°R 
i  =  2,  3,  .  .  . ,  200. 


192 


(4)  Calcu*.  tion  cl  6-r\ 


5(p)  =  (U3  r  Ay  *‘e3/5)/s4/5 


S  =  Sj  +  Peuet*e  r  Xj  initially.  At  x  +  Ax 


S  =  p  it  u  r  Ax 
x+Ax  e  e  e 


(5)  Calculation  of  v-derivatives  and  v-profile 


9v 

ay 

9M 

9y !: 


-e-) 


Ay 


u  -  u 
Hi+1  Pi-1 


2  Ay 


3u  j 


1 


(9%  ~2u„  “  6u2) 


3y|9  6Ay  '  3  "4 

1 


a2 
9  u 


*y 


ju 

dy 

d2I 

dy2 


9u 


91 

-5y 


A2 
o  u 

a  2 

ay 


6Ay‘ 


1 

6Ay 


2  (U4  "  3U3  -  9U2> 


i  =  2 


<9I3-2I4"6I2-V 


1 


2  6Ay‘ 


Ui+1  "  Ui-1 
2  Ay 


Ii+!  -  *1-1 


<I4-3I1-Sl9  +  5I1> 


2  '4 


2  Ay 


u.  +  u.  -  2u. 
l+l  l-l  i 


Ay 


193 


i  =  2,  3,  .  .  MIN  (j,  199) 


The  edge  of  the  boundary  layer  is  defined  by  j  =  i  at  which  the  following 
criteria  is  satisfied: 


Ay 


9u  9u  | 

9y  "  9y|e 


<  .00001 


<.  00005 

I  9y 


The  following  are  calculated  for  the  range  ^  <  i  <  j 


\  =■-  TT/(8u/ay) 
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PART  1. 

i 


2  (1  +  i 

Ui  '  r  dx  +  P  dx  / 


PART  2.  = 

l 


1 

P.u. 

l  i 


PART  3.  = - 

i  P.u. 
l  l 


dq 

+  (u:j 


1 

2Ay 


PART  6.  =  u. 

l  l 


2  v 

i-1  9v 

+  dy  j_j 


v.  =  - - PART  1.  +  u.  PART  2. 

l  du  2uj  I  li  l 

9v  Ay 
-  i  J 


- (PART  3 

(I  -  Uj  /2) 


-  u.  PART  2.)  -  PART  6. 

l  l  l 

dv  _  _2_  dv_ 

9y  j  "  Ay  v>  ~  Vr  "  ay  j.j 


195 


196 


u.  =  u. 

1  x+Ax  '  x+Ax 


■x+Ax  ^*x+Ax 

c.  Stagnation  Region  Calculation 


i  =  j  +  1,  j  +  2,  ....  200 


A  slight  modification  to  the  previous  discussed  procedure  is  used  to  obtain 
stagnation  point  profiles.  Input  profiles  are  corrected  by  integrating  the  u  and  I 
profiles,  but  x  (not  equal  to  zero)  is  not  increased  during  the  integration.  The 
profiles  are  assumed  correct  when 


‘  • 1  “list  for  a"  1 


This  convergence  criteria  is  obtained  with  the  velocity  similarity  stated  below. 


%  9u 
u  9x 


Experience  has  indicated  that  an  enthalpy  profile  convergence  criterion  is  not 
necessary. 

For  the  stagnation  region,  the  equations  for  the  velocity  and  enthalpy  profiles  are 
changed  to 


=  I  + 

i ,  S-AS 


9x 


AS 


9u 

u  +  — 

i,S-AS  9x 


x 

o 


x  +  AS 
o 


If  Ij  or  u^  become  negative  for  some  i,  AS  is  halved  and  the  case  is  restarted.  This 
may  be  repeated  a  maximum  of  two  times  for  any  combination  of  I,  or  Uj  <  0. 

Additional  input  quantities  are  AS  and  MITR. 

MITR  is  the  maximum  number  of  iterations  that  the  program  will  make  ir.  trying 
to  converge  on  a  solution  to  the  stagnation  problem. 
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AS  is  a  pseudo  length  similar  to  Ax.  The  program  also  uses  AS  as  a  test  to 
determine  whether  the  stagnation  point  calculation  is  to  be  performed  or  bypassed. 

For  a  flat  plate  case  AS  =  0. 

Calculations  for  a  hemisphere  may  be  performed  near  the  stagnation  point  and  also 
along  the  surface.  This  is  accomplished  by  the  input  of  AS  >  0  and  Ax  >  0. 


3.  INPUT-OUTPUT  DESCRIPTION 
a.  Data  Input  Preparation 


Data  are  input  through  the  medium  of  IBM  cards.  The  purpose  of  this  section  is 
to  define  the  required  inputs  and  the  form  they  take  on  the  cards. 

A  sample  data  sheet  may  be  found  at  the  end  of  this  section.  Use  of  this  type  of 
form  greatly  reduces  the  amount  of  effort  required  to  obtain  results  from  the  non¬ 
similar  boundary  layer  program.  The  data  may  be  keypunched  directly  from  the 
form. 

All  cards  except  cards  1  and  4  must  be  keypunched  with  decimal  points  and  power 
of  ten  indicators  as  follows:  if  only  a  decimal  point  is  needed,  the  number  may  be 
punched  anywhere  within  the  respective  ten  column  field;  however,  if  the  power  of  ten 
indicator  E  is  used,  the  exponent  must  be  right -adjusted  in  the  field. 

The  information  on  card  1  is  alphanumeric  and  may  be  keypunched  as  desired. 

The  data  on  card  4  are  all  integers  and  are  keypunched  as  fixed  point  numbers 
(i.  e. ,  no  decimal  or  power  of  ten  indicator  E  is  to  be  used).  Each  number  on  card  4 
must  be  right-adjusted  in  its  respective  five  column  field.  For  example,  if  JOT1  -  1 
then  a  one  must  be  keypunched  in  column  five. 


1)  Data  Card  Formats 


Card 

Columns 

F  ormat 

Description 

1 

1-48 

8A6 

Any  information  may  be  given  on  this  card.  It 
is  used  to  identify  the  case. 

2 

1-10 

7F10 

Ax,  the  increment  in  the  x  direction.  See  the 
discussion  on  numerical  stability  in  Appendix  C 
(Volume  I)  for  the  criteria  to  be  used  in 
choosing  a  value. 

2 

11-20 

7F10 

Ay,  the  increment  in  the  y  direction.  Should  be 

chosen  to  provide  at  least  ten  points  in  the 
boundary  layer  for  accurate  results. 
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Card 

Columns 

Format 

Description 

2 

21-30 

7F10 

T?max’  limiting  value  on  the  number  of 

calculations  in  the  y  direction. 

2 

31-40 

7F10 

Xq,  the  initial  value  of  x.  x0  >  0. 

2 

41-50 

7F10 

xf,  the  final  value  of  x. 

2 

51-60 

7F10 

<7,  the  turbulent  Prandtl  number. 

2 

61-70 

7F10 

CHECK,  case  dependence  criteria.  If  CHECK 
=  0  then  the  current  case  is  not  dependent  on 
previous  calculations.  If  CHECK  =  1  then  the 
current  case  is  dependent  on  the  previous  case. 

3 

1-10 

4F10 

AS.  Calculations  may  be  made  for  a  blui  t 
body.  This  involves  calculating  profiles  near 
the  stagnation  point  as  well  as  along  the  body. 
This  type  of  problem  is  indicated  by  AS  >  0 
and  MITR  >  0  (card  4).  For  nonstagnation 
point  calculations,  AS  =  0. 

3 

11-20 

4F10 

Turbulent  Prandtl  number  -  must  be  equal  to 
Prandtl  number  on  card  2. 

3 

21-30 

4F10 

Sj,  value  of  S  used  to  initialize 

S  =  S  +  p  u  r^  x 
lee  I 

3 

41-50 

4F10 

Xj,  value  of  x  used  to  initialize 

2 

S-S  +  pu  r  xT 
I  e  e  I 


The  data  on  card  4  control  the  number  of  values  to  be  input  to  each  of  the  suc¬ 
ceeding  tables.  The  values  must  be  greater  than  zero.  There  is  an  upper  limit  on 
the  number  of  values  that  may  be  input  to  each  table.  This  limit  is  given  with  the 
description  of  the  table.  The  limit  may  be  changed  by  altering  DIMENSION  statements 


in  the  program. 

4 

1-5 

1415 

JOT1, 

DJOT 

the  number  of  tabular  values  in  the 
and  XJOTT  tables.  0  <  JOTl  <  10. 

4 

6-10 

1415 

LIST, 
and  I 

the  number  of  tabular  values  in  the 
profile  tables.  2  <  LIST  <  50. 
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Card 


Columns 


Format 


Description 


4 


4 


4 


4 


4 


4 


5 


2 


6 


11-15  1415 


16-20  1415 


21-25  1415 


26-30  1415 


31-35  1415 


36-40  1^15 


NR,  the  number  of  tabular  values  in  the  three- 
dimensional  flow  parameter  table  (r),  the  table 
of  its  x-derivatives,  and  the  tables  of  their 
corresponding  values  of  x.  2  <  NR  <  50. 

NP,  the  number  of  tabular  values  in  the  pres¬ 
sure  table  and  the  table  of  corresponding 
values  of  x.  2  <  NP  <  50. 

NW,  the  number  of  tabular  values  in  the  wall 
enthalpy  table  and  the  table  of  corresponding 
values  of  x.  2  <  NW  <  50. 

NDU,  the  number  of  tabular  values  in  the 
(9u/9y)e  table  and  the  table  of  corresponding 
values  of  x.  2  <  NDU  <  50. 

NV1,  the  number  of  tabular  values  of  velocity 
normal  to  the  surface  and  the  table  of  corre¬ 
sponding  values  of  x.  2  <  NV1  <  50, 

MITR,  the  maximum  number  of  iterations 
allowed  to  calculate  the  profiles  for  the 
stagnation  region. 


Cards  5  and  6  are  used  to  specify  different  frequencies  of  output 
and  the  corresponding  ranges.  The  maximum  number  of  values  in 
these  tables  is  10. 


7F10  Table  DJOT,  assigns  value  to  first,  second, 

and  successive  printout  intervals. 

7F10  Table  XJOTT,  assigns  value  of  x  up  to  which 
each  interval  of  table  DJOT  is  used. 


2  In  the  case  of  tables  with  more  than  seven  input  values,  the  card  number  refers 
to  a  set  of  carls  (e.  g.  ,  if  ten  entries  are  made  in  the  DJOT  table,  then  card  5 
actually  refers  to  two  cards). 


Example:  Xq  =  0. 5  and  xf  =  1. 0.  Printout  is  desired  every  .  05 
(value  of  x)  until  x  =  .25,  every  .01  until  x  =  .3,  and 
every  .  1  until  x  =  xj.  The  required  input  would  be: 


Card  Column  Value 


i_4j 

; 

1-5 

3 

5 

1-10 

.05 

11-20 

.01 

21-30 

.1 

6 

1-10 

.25 

11-20 

.3 

21-30 

1.0 

(right-adjusted) 


Card  Columns  Format  Description 

7  -  7F10  Table  u-profile  =  f(y).  The  purpose  of  this 


table  is  to  provide  the  program  with  an  initial 
velocity  profile  at  x  =  xG.  Normally,  this 
table  must  be  input  and  contain  at  least  three 
values.  Exception  is  made  when  a  case  is 
dependent  on  a  previous  case  (see  Section  2) . 


This  table  contains  LIST  tabular  values. 

The  maximum  number  of  values  in  this  table 
is  50. 


7 

1-10 

7F10 

Velocity  at  wall. 

7 

11-20 

7  F 10 

Velocity  at  Ay  from  wall. 

7 

21-30 

7F10 

Velocity  at  2Ay  from  wall. 

8 

7F10 

Table  I-profile  =  f(y).  The  purpose  of  this 

table  is  to  provide  the  program  with  an  initial 
total  enthalpy  profile  at  x  =  xQ.  Normally,  this 
table  must  be  input  and  contain  at  least  three 
values.  Exception  is  made  when  a  case  is 
dependent  on  a  previous  case  (see  Section  2). 

This  table  must  contain  the  same  number  of 
values  (LIST)  as  the  u-profile  table.  The 
maximum  number  of  values  in  this  table  is  50. 
Entries  are  assumed  to  be  consistent  w  ith  Ay. 

8  1-10  7F10  Total  enthalpy  at  the  wall. 
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Card 

Columns 

Format 

Description 

8 

11-20 

7F10 

Total  enthalpy  at  Ay  from  the  wall. 

8 

21-30 

7F10 

Total  enthalpy  at  2Ay  from  the  wall. 

9 

7F10 

Table  r  =  f(x).  This  table  contains  NR  tabular 
values  of  r,  the  three-dimensional  flow  param¬ 
eter  (streamline  divergence  due  to  body 
geometry).  The  maximum  number  of  values  in 
this  table  is  50. 

10 

7F10 

Table  (dr/dx)  =  f(x).  This  table  contains  NR 
tabular  values  of  (3r/3x).  The  maximum 
number  of  values  in  this  table  is  50. 

11 

7F10 

Table  x.  This  table  contains  NR  tabular  values 
of  x  corresponding  to  entries  made  in  the  r 
and  (3r/3x)  tables  (cards  9  and  10). 

12 

7F10 

Table  pressure  =  f(x).  This  table  contains  NP 
tabular  values  of  pressure.  This  table  is 
numerically  differentiated  for  (3P/3x).  The 
first  value  of  (3P/3x)  is  set  to  zero.  The 
maximum  number  of  values  in  this  table  is  50. 

13 

7F10 

Table  x.  This  table  contains  NP  tabular  values 
of  x  corresponding  to  the  entries  made  in  the 

P  table  (card  12). 

14 

7F10 

Table  wall  enthalpy  =  f(x).  This  table  <  ntains 
NW  tabular  values  of  the  total  enthalpy  c  .he 
wall.  This  table  is  used  for  nonisothermal 
wall  problems.  The  maximum  number  of 
values  in  this  table  is  50. 

15 

7F10 

Table  x.  This  table  contains  NW  tabular  values 
of  x  corresponding  to  the  entries  made  in  the 
wall  enthalpy  table  (card  14). 

16 

7F10 

Table  (3u/  3y)e  =  f(x).  This  table  contains  NDU 

tabular  values  of  (3u/3x)  evaluated  at  the  edge 
of  the  boundary  layer.  For  all  problems  not 


involving  vorticity  (3u/3x)  =  0.  The  maximum 
number  of  values  in  this  table  is  50. 
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Card 


Columns 


Format 


Description 


17  -  7F19  Table  x.  This  table  contains  NDl'  tabular 

values  of  x  corresponding  to  the  entries  made 
in  the  (du/9y)e  table  (eard  Hi). 

18  -  7F10  Table  Vw  =  f(x).  This  table  contains  NV1  tab¬ 

ular  values  of  the  \elocity  normal  to  the  wall. 
This  table  is  used  for  mass  injection  or  leakage 
problems.  The  maximum  number  of  values  in 
this  table  is  50. 


19  -  7F10  Table  x.  This  table  contains  NY1  tabular 

values  of  x  corresponding  to  the  entries  made 
in  the  Vw  table  (card  18). 

20  NOTE:  In  all  cases  where  the  values  in  the  above  tables  are  zero, 

a  zero  must  be  entered  on  the  input  sheet. 

Example:  the  normal  velocity  at  the  wall  is  zero: 

Vw  table  (eard  14)  0.0  0.0  0.0 

x  table  (card  15)  0.0  1.0  2.0 

2)  Case  Dependence 

The  NSBL  program  has  a  criterion  which  allows  a  case  to  be  dependent  on  the 
suecess  or  failure  of  the  previous  ease.  This  involves  the  input  CHECK  and  LIST 
and  the  test  for  ERROR. 

If:  Then: 

LIST  CHECK  ERROR 

0  0  —  Not  dependent  on  the  previous  case,  program 

continues. 

0  1  0  Dependent  on  the  previous  case  which  worked, 

program  continues,  using  final  profile  from 
previous  case. 

0  1  1  Dependent  on  the  previous  case  which  failed,  so 

present  case  stops  and  program  continues  to 
next  case. 

0  -  0  Dependent  on  the  previous  case  which  worked, 

program  continues,  using  final  profile  from 
previous  case. 

Dependent  on  the  previous  case  which  failed, 
so  present  case  stops  and  program  continues 
to  next  case. 
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0 


1 


b.  Output  Description 

The  printed  output  consists  of  the  input  data  and  at  each  printout  interval  the 
following  are  printed  as  a  function  of  y: 

y  u  q  3v/3y  v  3u/3x 

T  t  I  17  3I/3x  e 

The  following  are  printed  as  a  function  of  the  x  at  the  end  of  each  printout  interval: 

Tw,  L 
tw,T 
a* 
e 

\ 

3P/ 3x 

QQW 

S 
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4.  PROGRAMMING  INFORMATION 
a.  Program  Flow  Chart 


*J 

INPUT 

n 

_ I _ L 

PRINT  INPUT  j 

dTj, 


LOOK  UP 
_  9P  9r  9u  | 

P  9x'r'  9x*  9vle'vw 
AS  FUNCTIONS  OF  x 

I 

CALCULATE 

an  du  91  a2 U  92l  arL  9tT  ,  T* 

‘  9V  '  9V '  8y2  '  ay2  ’  9v  '  9y  '  5- '  T ' 

FOR  ALL  POINTS  IN  THE  PROFILE  IT)  <  T)^) 

I 

CALCULATE 

T.  .  tt.  t<  e,|^,  AO,  O’,  PART  1,  PART  2.  PAST  3 
L  1  oV 

PART  6.  w,  FOR  ALL  POINTS  IN  PROFILE 

9v  9x 


CALCULATE  qw,  Tw,  <5#,  9 


i  PRINT  OUTPUT  1 

_ : 1  I  - 


INCREASE  x  lYAx  I 

- - 


LOOK  U? 

p  9P  9r  j  9u 

'  9x  '  ,r  c  x'  w'  9y 
AS  FUNCTION  OF  x+i 

-V1 

e 

^x 

_ i _ 

CALCULATE  u  AND  1  PROFILE 
p,  T,  p  AT  x  +Ax 

2 


I 


b 


NSBL  Program  Listing 
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